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Abstract 

This  thesis  developed  a  finite  element  solution  of  the 
Schrodinger  wave  equation.  This  technique  is  used  by  a 
computer  program  to  calculate  the  energy  levels  and  wave 
functions  of  a  diatomic  molecule  for  a  particular  poten¬ 
tial  energy  model.  The  potential  energy  model  is  a  function 
of  a  set  of  parameters  which  a  non-linear  minimization  rou¬ 
tine  varies  before  solving  the  wave  equation.  This  is  done 
in  an  iterative  manner  until  the  calculated  energy  levels 
agree  in  a  least  squares  sense  with  the  observed  energy 
levels.  Then  the  transition  probabilities  ( Franck-Condon 
factors)  between  the  wave  functions  are  calculated  by 
another  program  developed  for  this  thesis.  Finally,  two 
programs  were  written  to  determine  the  energy  levels 
observed  in  spectroscopic  data.  One  uses  Dunham  coefficients 
and  the  Dunham  equation  while  the  second  uses  a  least  square 
fit  to  the  data  directly. 

The  four  programs  were  tested  and  appear  to  work  correct¬ 
ly.  The  numeric  solutions  were  compared  with  the  analytic 
solutions  of  the  single  harmonic  oscillator.  The  lowest 
2‘j  energy  levels  agreed  to  within  0 . 005#  accuracy  while 
their-  wave  functions  appear  to  agree  to  within  0.40% 


COMPUTER  MODELING  OF  VIBRATIONAL  ENERGY  LEVELS 

OF  POTENTIAL  LASER  CANDIDATES  (DIATOMIC  MOLECULES) 

1 .  Introduction 

Back. 9  round 

The  evaluation  of  the  lasing  potential  of  diatomic 
molecules  is  simplified  by  accurate  knowledge  of  the 
molecule's  energy  as  a  function  of  internuclear  distance. 
This  knowledge  is  used  to  generate  wave  functions  from  the 
Sohrodinger  wave  equation  which  describe  the  molecule  in  a 
particular  quantum  state.  These  wave  functions  are  n. 

turn  to  calculate  the  probability  that  the  molecule  W’  11 
change  from  one  state  to  the  other.  This  probability,  or 
Fruuck-Ccndon  factor,  eases  the  correlation  of  spectroscopic 
data . 

AFIT  began  an  effort  to  develop  a  set  of  computer 
routines  capable  of  calculating  these  Franck-Condon  factors 
for  diatomic  molecules  in  1982.  The  central  program  of  this 
set  was  acquired  from  Dr.  C.R.  Vidal  (Max  Planck  Institut  fur 
Extruterrestriche  Pt.ysik)  (26).  This  program  uses  the  semi- 
classical  Kydberg-Klein-Rees  (RKR)  procedure  and  an  Inverse 
Pe r t urba tiona 1  Analysis  (IPA)  to  generate  the  molecule's 
potential  energy  curve.  Capt.  L.L.  Rutger  modified  Vidal's 
program  to  run  on  the  CDC  CYBER  computer  system  in  his  thesis 
effort  (March  83).  Rutger  also  wrote  a  program  to  generate 


t  he 


set.  of  molecular  constants  used  as  input  to  the  RKk-IPA 

(20).  Then,  in  another  thesis  project  (December  83), 
Capt.  J.J.  Pow  compl e t.ed  the  program  sl  t  by  creating  programs 
to  plot  the  energy  curves  and  calculate  the  Fr a nck-Condon 
laciors  (It). 

The  RKR-iPA  program  set  yielded  results  that  agreed 
very  well  with  similar  work  in  the  literature  (5).  However, 
interest  uoes  exist  in  finding  a  more  efficient  program  set 
based  on  some  oilier  technique*  than  RKR.  This  technique 
should  avoid  the  following  shortcomings  of  the  RKR  program 
set.  First,  the  RKR- 1 PA  programs  require  considerable 
com£  ut.ei  resources  and  arc-  therefore  expensive  to  use.  Tnese 
programs  can  not  run  on  just  any  minicomputer  that  may  be 
available.  Tney  require  full  mainframe  support.  Second, 
researchers  have  observed  anomalous  behavior  of  the  RKR 
potential  (If),  ( i  0 )  ,  (  2  3),  (2  5),  (29  ).  The  curve  sometimes 
be itd a  over,  or  decreases  rapidly  for  higher  energy  states. 
These  effects  are  usually  attributed  to  an  incorrect  set  of 
molecular  constants.  Wells,  Smith,  and  Zare  concluded  that 
the  RKR  method  is  vt ry  sensitive  to  errors  or  inconsistencies 
in  the  experimental  data  (27).  This  has  sparked  interest  in 
u  t. echniquo  tiidt  uses  the  experimental  data  directly  to  find 
the  potential  that  yields  the  best  fit  to  the  measured  energy 

L  1  i  1 1  i  1 » C  t  ;  Li  • 

Previous  work  to  find  a  new  technique  or  better  version 
ui  RKR  have  not  made  any  significant  improvements  over  RKR, 


Now  the  vector  can  be  back-transformed  into  vector  x, 
which  contains  the  Values  of  the  energy  levels,  by  solving  LT 
±  -  X*  L  pp  -  0  since  L  and  L  have  the  same  diagonal.  This 
weans  that  Xp  is  urbritary  since  0  xp  =  Yp  where  yp  is  also 
arbitrary.  The  vector  element  Xp  is  chosen  to  be  zero  since 
it  at  loots  the  values  of  other  vector  elements  (i  <  p ) . 
The  remaining  elements  ot  ^  are  :  elated  to  x  by 


T 

L.  -  x. 

li  l 


P-1  m 
r  a  L*  x 
k=i+llK  K 


(i  1 i 2, ...» (p  1)) 


since  Lt  is  upper  triangular.  Inverting  and  remembering  I^jj. 
-  Lkl  Eg  (14)  becomes: 


P_1 

—  Yy  “  ^  (  i "( P“1 )»( P“2 1 )  (15) 

1  k  i+l  K 


The  value  of  :<p_^  is  computed  first  since  it  is  needed  to 
compute  the  value  of  x  2 • 

The  vector  x  now  contains  the  energy  levels  that  best 
fit  the  spectroscopic  data.  However,  the  highest  energy 
level  '■  has  a  value  of  zero  and  all  other  levels  are 
negative.  The  energy  levels  may  now  be  shifted  so  that  t  ^  is 
either  zero,  or  the  value  the  Dunham  equation  (Eg  (1)) 
yie  ids . 

This  least  squares  method  yields  values  for  only  the 
energy  levels  represented  in  the  data.  Missing  values  may  be 
found  from  the  Dunham  equation  (Eq  (l)). 


( j-1) 

- 


(j-l)o  A 

=  (A,,  -  £  L?v )2 


where  L ^ j  =  0  if  j  >  i  since  L  is  lower  triangular.  Also, 
the  last  diagonal  element  Lpp  is  zero  since  A  is  singular. 

The  solution  of  L  y  =  b  is  found  by  inverting 

( i-1 ) 

bi  =  “lA  +  Liiyi 

K  — 1  /  i  1  v 


to  9et 


(i-D 

"  bi  '  ALikyk 


This  method  yield.-,  only  p  -  1  values  for  \/  since  Lpp  =  0  and 
i-;q  (ID  call  not  be  solved  for  y  This  means  that  yp  is 

arbitrary  and  can  be  chosen  to  be  zero.  However,  the 


cond i t l on 


(p-1 ) 

»P  = 


must  be  met  it  the  problem  is.  to  be  consistent 


level  difference  information.  The  solution  is  unique  only  in 
terms  of  differences  not  in  the  value  acquires. 

Since  A  is  positive  semi-definite  (28:28-30),  the 
Cholesky  decomposition  (28:229)  (4:Sec  8  1-12)  can  be  used  to 
state  the  problem  in  a  solvable  form.  First,  a  lower 
triangular  matrix  L  is  found  so  that  A  =  L  LT  (L ^  is  the 
transpose  of  L).  Then  the  problem  becomes  L  ^  =  b  where  = 
Ltx. 

The  decomposition  of  A  occurs  as  follows.  If  A  =  L  l/1  , 


then  element  A^j  of  A  is: 


min(i, j)m 

t  t 


=  L 
k=l 


LikLicj 


min( i, j) 

=  k^LikLjk 


(8) 


where  LTkj  =  Ljk  since  LT  is  the  transpose  of  L.  Three  cases 
arise  when  evaluating  Eg  (8): 


P  (j-l)p 


Cj 


(j"l) 

LijLjj  +  k=1LikLjk 


(9) 


1<J 


Aij 


(i-1) 

LijLii  +  k^LikLjk 


Inverting  Eqs  (9)  yields 


12 


0 


(6) 


=  J/lk'hk'W  _jl1Wkj<1kj"€k','Gj> 

A  final  rearrangement  of  terms  yields: 


t-J1<wik+wki)  -J1(wkrwjk)<:j  =J1(wkjikj-wjkijk) 


k. 


(7) 


The  first  term  of  Eq  (7)  is  the  sum  of  tne  weighting  factors 
for  all  transitions  ending  at  level  £  ^  plus  the  sum  of  the 
weighting  factors  for  all  transitions  beginning  at  level  £  ^ 
times  the  value  of  The  second  term  is  the  sum  over  j  of 
the  weighting  factor  for  the  transitions  between  levels  £  ^ 
and  €_j  times  the  value  of  Finally,  the  third  term  of  Eq 
(7)  is  the  sum  over  j  of  the  weighted  transition  from  level 
1  k  to  level  minus  the  weighted  transition  from  level  £  j  to 
level  0^.  Only  observed  transitions  are  used  in  Eq  (7) . 

This  relationship  can  be  more  conveniently  solved  by 
rewriting  Eq  (7)  in  matrix  form  as  shown  in  Figure  11-3. 
This  is  the  common  linear  problem  Ax  =  b  where  A  is  a  p  by  p 
square  matrix,  and  the  vectors  x  and  b  are  p  by  1  matricies. 


The  diagonal  elements  of  A  (A 


—  '"inn 


m 


=  n)  are  the  sum  of  (w^  + 


w 


K  1 


i  =  1  to  p. 


The  off  diagonal  elements  (A 


Bin 


m 


/  n)  are 


~^wmn  +  wnm*' 


The  problem  is  complicated  by  tire  fact  that  A  is 
singular,  possessing  no  inverse  A” ^  such  that  x  =  A--'-  b.  A 
is  singular  since  the  spectroscopic  data  contain  only  energy 


10 


inconsistencies  in  the  additional  data.  The  following 
analysis  finds  the  set  of  energy  level  values  that  best 
resolves  the  inconsistencies  of  the  entire  set  of  data. 

The  weighted  linear  least  squares  technique  finds  a  set 
of  energy  levels  with  the  smallest  sum  5  for  the  p  observed 
levels . 


S 


P  P 

i  £  £ 

i=l j=l 


w-  •  ( 1  • 

iy  ij  i  j 


(2) 


The  weighting  factor  w^j  allows  more  accurately  determined 
data  to  make  a  larger  contribution  to  the  sum  S.  The  minimum 
of  the  sum  S  occurs  when  its  derivatives  with  respect  to  the 
energy  levels  of  interest  are  zero: 


The  derivative  is 


is  a  o 


(3) 


0  = 


P 

£ 


P 

£ 


i=l j=l 


wij(1ij-6i+5j)<6jk-6ik) 


(4) 


since 


d  € 

aT 


m  = 


ran 


n 


1  if  m=n 
0  if  m^n 


(5) 


Using  this  property  (Eq  (5))  of  the  Kronecker  delta  6mn,  Eq 
(4)  becomes: 
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coefficients  are  often  published  in  place  of  the  data  they 
represent.  The  Dunham  coefficient  Y^j  is  a  coefficient  of  an 
infinite  power  series  of  vibrational  and  rotational  quantum 
numbers  (v  and  J).  Given  enough  terms,  this  series 
accurately  represents  the  energy  T(v,J)  of  the  quantum  state 
represented  by  v  and  J  (6:725).  The  Dunham  equation  is: 


T(y  j)  =  Xl_Yi.(v+^)1JJ(J-KL)J 


Pow  describes  a  good  method  of  determining  the  Dunham 
coefficients  from  spectroscopic  data  (16:4-8). 

Alternatively,  the  researcher  may  use  the  spectroscopic 
data  directly  to  determine  the  observed  energy  levels  by  a 
weighted  least  squares  fit. 

First,  a  set  of  transition  lines  l^j  are  observed  and 
assigned.  The  assignment  identifies  the  transition's  initial 
and  final  electronic  and  vibrational  states.  Each  transition 
is  from  an  initial  energy  level  £ j  to  a  final  level  £j. 
These  levels  may  or  may  not  belong  to  the  same  electronic 
state  as  shown  in  Figure  II-l.  For  example,  line  I73 
represents  a  transition  from  electronic  state  B,  v  =  2  to 
state  A,  v  =  2.  Also,  line  196  represents  the  transition 
from  state  C,  v  =  2  to  state  C,  v  =  0. 

Then,  the  least  squares  technique  is  used  to  resolve 
inconsistencies  between  the  data.  A  transition  line 


involving  one  energy  level  may  not  agree  with  other 


Literature  Background  and  Theor 


1 

This  section  develops  concepts  and  theory  required  to 
understand  the  programs  presented  later.  The  goal  of  this 
effort  is  the  calculation  of  the  transition  probability 
(Franck-Condon  factor)  between  two  energy  states  of  a 
diatomic  molecule.  The  wave  functions  describing  these 
states  are  used  to  calculate  the  Fr anck-Condon  factor,  and 
are  found  by  solving  the  Schrodinger  wave  equation.  However, 
knowledge  of  the  molecule's  kinetic  and  potential  energy  is 
required  before  solving  for  the  wave  functions.  Therefore,  a 
discussion  is  first  presented  on  determining  the  appropriate 
electronic  and  vibrational  energy  levels  of  the  system. 
Then,  two  potential  energy  models  are  presented.  This 
section  is  followed  by  a  detailed  outline  of  a  finite  element 
solution  of  the  wave  equation.  Finally,  the  calculation  of 
Franck-Condon  factors  is  discussed  (21). 

determination  of  Energy  Levels 

There  are  two  starting  points  for  determining  the  energy 
levels  of  a  diatomic  molecular  system.  The  researcher  may 
either  use  Dunham  coefficients  to  calculate  the  energy 
levels,  or  they  may  be  derived  directly  from  spectral  data. 

Dunham  coefficients  are  combinations  of  molecular 


constants  derived  from  spectroscopic  data. 


These 


I  iF'Ji»g  u  p  in 


not  be  important,  as  long  as  it  allows  the  flexibility  to 
modify  the  repulsive  and  attractive  branches  independently. 

This  program  was  validated  against  the  analytic 
solutions  of  the  simple  harmonic  oscillator.  Two  programs 
have  been  written  to  aid  tlie  researcher  in  defining  the 
energy  levels  used  as  input.  One  program  uses  a  least 
squares  technique  to  find  the  set  of  energy  levels  that  best 
fit  a  set  of  spectroscopic,  transition  data.  The  second  uses 
Dunham  coefficients  to  calculate  approximate  values  for 
missing  energy  levels.  The  last  program  written  computes  a 
25  by  25  (v*  -  0  to  24;  v"  =  0  to  24)  table  of  Franck-Condon 
factors  between  two  electronic  states. 
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Choose  a  set  of  vaules 
for  the  potential  energy 
parameters  a,  ... 


Solve  the  wave 
equation: 


Ecalc  = 


Compute-. 

S  =  ilw(Eobs  -  Ecalc>2 


Fig.  1-i. 


Minimization  of  the  Least  Squares  Sum 


or  are  still  computationally  expensive  (8),  (12),  (24). 


Objective 

The  objective  of  this  work  is  to  develop  and  validate  a 
program  set  based  on  a  quantum  mechanical  approach.  These 
programs  should  be  transportable  to  minicomputers  such  as  a 
VAX  11/780,  HP  1000,  Harris  800  and  be  inexpensive  to  use. 
Also,  they  should  use  observed  energy  levels  as  input,  not  a 
set  of  molecular  constants  derived  from  these  observed 
levels . 

Approach 

A  program  has  been  written  which  solves  the  wave 
equation  using  the  finite  element  method.  This  program  uses 
a  potential  energy  function  defined  in  terms  of  parameters 
the  program  can  change.  A  non-linear  minimization  routine  is 
used  to  find  the  set  of  parameters  which  characterize  the 
potential  that  best  fits  the  observed  energy  levels  in  a 
least  squares  sense.  This  approach  is  illustrated  in  Figure 
1-1.  The  minimization  routine  selects  a  set  of  parameters, 
solves  the  wave  equation,  and  computes  a  sum  of  the 
difference  between  observed  and  calculated  energy  levels 
squared.  This  process  is  repeated  until  the  set  of 
parameters  is  found  that  yields  the  smallest  energy 
difference  sum. 

This  approach  allows  the  researcher  to  fit  a  potential 
energy  model  to  his  data.  The  specific  model  chosen  should 


The  molecular  constant  £  is  defined  as  (10:101): 

@  =  1-217?, 

where  pi  is  the  reduced  mass  of  the  system  and  (Jfc  is  a 
vibrational  frequency.  Tel  1 inghuisen  and  Henderson  observed 
that  the  Morse  function  seems  to  accurately  describe  the 
repulsive  branch  ot  the  potential  (24). 

Another  potential  energy  model  is  the  Lennard- Jones 
potential : 


This  potential  is  considered  since  it  is  the  sum  of  a 
repulsive  and  attractive  potential.  Therefore,  each  branch 
may  bo  modified  independently  of  the  other  branch. 

Whichever  model  is  chosen,  the  potential  is  described  in 
terms  of  parameters  than  can  be  altered  before  solving  the 
wave  equation.  For  example,  the  powers  12  and  6  of  Eq  (IB) 
might  be  replaced  by  a  and  pi  .  Then  a  and  pi  are  each  allowed 
to  have  a  certain  value,  say  a  =  12  and  (3  =  6  or  a  =  11.39 
arid  0  -  5.06.  The  method  then  is  to  vary  the  potential 
parameters  and  solve  the  wave  equation  in  an  iterative 
fashion.  Given  enough  iterations  the  set  of  parameter  values 
will  be  found  for  the  potential  energy  model  making  the  best 
fit  to  the  observed  energy  levels. 


Care  must  be  exercised  when  the  potential  energy 
function  is  parameterized  to  ensure  that  the  model  is  still 
an  accurate  description  of  the  molecule.  For  example,  if  the 
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Lennard-Jones  potential  were  parameterized  as 


(19) 


one  would  find  that  the  function's  minimum  no  longer  occurs 
at  re  if  a=  11.39  and  $=  5.06.  However,  if  the  potential 
were  parameterized  as 
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then  the  function  behaves  well  for  all  a  and  (3.  The  function 
described  by  Eg  (20)  is  the  Mie  potential  (14:311).  Moelwyn- 
Hughes  points  out  that  both  the  Morse  and  Lennard-Jones 
functions  are  special  forms  of  the  Mie  function  (14:311-315). 


Numerical  Solution  of  the  Wave  Equation 

The  wave  functions  used  to  calculate  Franck-Condon 
factors  are  derived  by  solving  the  Schrodinger  wave  equation. 
The  procedure  uses  the  finite  element  method  to  solve  for  the 
energy  eigenvalues  of  a  particular  set  of  potential  energy 
parameters.  The  calculated  energy  eigenvalues  are  compared 
with  observed  energy  levels  in  a  weighted  least  squares 
sense.  Then,  a  new  set  of  potential  energy  parameters  is 
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generated,  the  wave  equation  solved,  and  energy  levels 


compared.  This  process  is  repeated  until  the  parameter  set 
that  best  fits  the  measured  energy  levels  is  found.  Then  the 
eigenvectors  are  computed  and  normalized. 

The  expected  values  of  the  energy  operator  are: 

E  =  extlj'dx^wi  (21) 
[j'dxy  i|f  J 

where  ext  (  1  means  the  extremum  of  the  backeted  quantity. 
The  extrema  of  a  function  are  its  minima,  maxima,  and  saddle 
points.  Assuming  spherical  symmetry,  Eq  (21)  becomes 


E  =  ext< 
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where  the  Laplacian  operator  is: 

V2  =  L_  d_  r2  d_ 
2 

r  dr  dr 


(22) 


(23) 


The  wave  function  (r)  is  replaced  by  a  function  U|rj  so 
that  : 
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(See  French  (7:199-201)  for  another  discussion  on  this 


substitution) .  Then  Eq  (22)  becomes: 
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the  kinetic  energy  term  of  Eq  (25)  is  rewritten  as: 


(27) 


Therefore  Eq  (25)  is  now 


E  =  ext<0 


JodrU<r>(1^V>)V> 


(28) 


When  integrated  by  parts,  Eq  (28)  becomes 
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J‘drU^ 


The  finite  element  method  can  now  be  used  to  solve  Eq 
(29).  First  a  uniform  grid  is  overlayed  on  the  wave  function 
as  indicated  in  Figure  II-5.  Each  grid  point  is  a  node. 
The  natural  coordinate  system  (Fig  11-6)  (3:88)  simplifies 
the  problem.  Each  position  r  is  described  in  terms  of  the 
local  grid  boundaries  r  ^  and  ri  +  1,  and  the  natural 
coordinates  1^  and  I2  as: 


r  =  q  +  l2h  =  ri+1-  qh 


where  h  is  the  grid  element  size 


h  =  ri+r  ri 


and  1^  +  1 2  -  1  • 

The  approximate  solutions  of  Eq  (29)  found  by  the  finite 
element  method  will  converge  to  the  true  solution  as  the  grid 
element  size  is  made  smaller.  However,  the  interpolating 
polynomial  chosen  to  approximate  U(rj  must  satisfy  the 
requirements  of  completeness  and  compatibility  presented  by 
hao  (17:114-118).  These  requirements  are  met  if  a  basis  set 
of  terms  cubic  in  1  ^  and  I2  are  used  to  approximate  Uir). 


Then  the  requirement  is  that  the  approximating  function  and 
its  first  derivatives  are  continuous  at  each  node  r^.  The 
function  U(r)  is  approximated  by: 


U(r)  =U(lrl2)  rfrsrm 

U(l1,l2)=U0fl(ll,l2)+t)0f2(l1,l2)+Ulf3(l1,l2)+Ulf'»(l1.l2) 
where  UQ  =  value  of  at 

I 

=  slope  of  at 

=  value  of  at 

I 

=  slope  of  U,j  at  ri+1 


(32) 


Then  the  boundary  conditions  at  the  left  edge  of  the  grid 
element  (  1^  -  1,  I2  =  0)  are: 


fl(l,0)  =  1 
f2(l,0)  =  0 
f3(l,o)  =  0 
f4(l,0)  =  0 


0 


f4(l,0)  0 


(33) 


The  boundary  conditions  at  the  right  edge  (  1 1  =  0 , 


12  =  1) 
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The  function  Ujrj  should  be  zero  and  have  zero  slope  at 
r  =  0  and  r  =■ 00  .  This  condition  is  not  enforced,  but  should 
be  a  result  of  this  analysis  as  the  function  gets  far  from 
the  potential. 

Tile  basis  set  for  cubics  in  1^  and  I2  is  £  1  ^ ;  1j_212; 
ljlo^;  l-,^j.  The  first  derivatives  of  the  basis  set  are  [- 
M1'//\r,  (-211l2  +  li2)/h;  (- 1 22  +  2  1  1  1 2) /h;  3122/h]  since 


and 


dy  =  dy  dli  +  dy  dig 
dr  dl  1  dr  dl  2  dr 


dli  =  -1 
dr  h 

(36) 

dls  =  1 

dr  h 


so  that  Eq  (35)  becomes 


dy  =  lfdx  ~  d^  \  (37) 

dr  h\dl2  dlj 


The  combination  of  basis  functions  that  satisfy  the  boundary 
conditions  (Eqs  (33)  and  (34))  is: 
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=  -(61^  )/h 


hllX2 


1&h1l 


-hlA 


=  (l2-2l1l2)/h 
=  (6l1l2)/h 


(lj-2l1l2)/h 


«.-■  function  u(r)  anc*  its  first  derivative  are  now 
P r ox i mated  by: 
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Ui.ip  =  U0(  •  3lfl2)  +U0(hlfl2)*U1 


+U1(-hl11|) 
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=  U0i^(3u0+U0h)iji2+(3U1-U1h)i1i2+U1i3 


=  u0^  6ll12yUp(-211l2tl^)-m1^0ll12 


^(-21^2-^lp 

=  Uq12+( -3u0/h-UQ+3U1/h-U^ ) 21112+U^li 


The  integrals  of  Eg  (29)  may  now  be  rewritten  as 


I  =  — 
1  2 


/di1di2[u0lj+(3u0+u’h)l2i2^(3u1-u^h)i1i^ 


(42) 


(43) 


h  3  2  2  3  2 

i3-Xdi1di2Lu0i1+(3u0+UQh)i1i2+(3u1-u^h)i1i2+u1i2] 


for  each  grid  element  where 


I«  +1 0~S  dr 
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These  integrals  are  solvable  by  the  simple  relation 
(3:312) : 


rdi  di  ipiq 

Jndlldi21li2  Tp+q+1)! 


(45) 


Since  U0,  UQ',  U  2  ,  and  U1'  are  constants  with  respect  to  the 
coordinates  1  3  and  1 2  /  the  result  is  the  integral  of  a 
polynomial  in  terms  of  1 1  and  12-  Then,  if  each  integral  is 
separated  into  a  sum  of  integrals  and  constants  are  taken  out 
of  the  integral,  Eq  (45)  solves  each  term  of  the  polynomial. 
The  integral  I2  (Eg  (42))  is  complicated  at  this  point  by  the 
potential  energy  term  V ( 1 x , 1 2 )  and  will  be  discussed  next. 
When  solved,  1-^  (Eq  (41))  and  I3  (Eq  (43))  result  in  a 
polynomial  in  the  boundary  values  UD,  U0*,  U^,  and  U^'  of  the 
grid  element.  For  example,  I3  becomes 
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16?2hU2  +  264h2UnU*  +  648hUnlL  +  . . . 
?!  7!  u  ?!  u  1 


This  expression  is  more  convenient  when  written  in  matrix 
i  o  r  hi  as: 
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Similarly,  the  integral  1 ^  becomes 
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The  integral  over  the  entire  grid  is  then  reduced  to  a  matrix 
problem  composed  of  one  matrix  I  j  for  each  grid  element  i. 
lor  example,  the  overlap  integral  becomes: 


^III  ^  dt'U/  )  —  .z  S  z_ 
a  ' 


where  a  and  b  are  the  first  and  last  nodes  of  the  grid.  Then 
zJS  L  is  constructed  as  illustrated  in  Figure  II-7  where  the 
j  *  *'  sub-mat  ri cies  of  z  and  S  are  shown  in  Figure  11-8.  bach 


The  value  and  slope  of  U,  >  at  the 


S- 

— 1 


I872h  2o4h2  X  X 

71  71 


2o4h2  X  XX 

7! 


X  XXX 


X  XXX 


Fig.  1 1-8.  The  Sub-Matricies  of  zTS  z 


sub-matrix  corresponds  to  one  grid  element.  The  sub- 
matricies  zl  and  S1  overlap  with  sub-matr icies  z_1-^ ,  z^  +  -*-  and 
F  +  1  since  the  (i-l)th  end  i1*1  grid  elements  share  the 
same  grid  boundaries,  us  do  the  (i+l)t^  and  i^1  elements: 
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(50) 


When  element  F|un  of  S  or  zm  of  z  correspond  to  the  overlap  of 
sub-mat  r  icies  b1,  +  l  or  rz}- ,  z^  +  ^  the  value  of  Smn  or  zm  is 

the  sum  of  the  corresponding  values  of  the  sub-matricies. 
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The  integral  (Eg (42))  is  solved  in  the  saii.e  manner; 
however,  the  potential  energy  function  V(rj  must  first  be 
approxin.ut  ed  by  VUj,^).  The  potential  V^r)  is  written  as  a 
combination  of  the  same  cubic  basis  set  used  for  the  function 


U|rj.  This  keeps  the  same  level  of  accuracy  between  the 
approximating  1  unctions  11(1^, 12)  and  V(l^,l2).  Specifically, 


V|t.j  becomes: 


V(r)-V(l1,l2) 


r  d 


ri-tl 


(SI) 


V(i1,i2)  vofl(l1.l2)+vof2(l1,l2)+vlf3(i1,i2)+Vlf4(l1,l2) 

where  the  ‘s  are  the  same  functions  used  earlier  (Eq  (38)). 
Then  the  integral  I2  (Eg  (42))  for  the  it^  grid  element 
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r2tp(r)V(r) 

=Xdl1dl2[U0lJ-t(3U0+UQh)l^l2+(3U1-U^h)i1l2+u1l3]2[v0lJ  (52) 
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(53) 


After  the  same  manipulations  described  before,  the  l2a's  are 
evaluated  as: 
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File  COMP 


Ml 

GE  SPOCK 

JE\2150  '  BBBA 

1  r-'.'EA  CO  VLCN2  SPOCK 

ED  SPOCK 

TA 

E.3 

C  ,  6  -  e ,  L 2 
*JE,  120,  1  MED 
1  NED  UP 
MO  *A-GH 
IK  0,o 

*3i  ,  ::?3,  1  NEXT 

■NEXT  *SAUF77.  LPUW&O  LI 

*JE  ! F77ERR 

*PR  compiled  ok... 

*JU  'ELOLD 
IF77ERR  *PR  OoDpt,... 

‘  ELOL D  EL  L2 
*  J  E , 312,  ! EL2Z2 
*JE,2170, ! NEXT2 
■*JU  !  NEXT 2 
1 EL22Z  EL  222 
KJE, 2170, ' NEXT 2 
1 MEXT2  JS  SPOCK 
110  S-A-OFF 
EL  SPOCK 
AS  a  -  * 

*PK  done  arid  did. 

ME 


File  VLCN2 


4>PR  vijl  1  Cart  I  Zat  i  or. .  .  . 

VU.  b 

NAME  222 

L  I  ,  2039EFPH*  BINDER  ,  1  OOOAFI  T  *  Il'ISL  1  B  ,  *SAUL77  , *LIBERY 
BE 


ig.  Ill—] .  Fortran  77  Compile  and  Link  Macro 
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The  job  control  statements  required  to  compile  and  link 
these  programs  is  given  in  Figure  III-l.  These  statements 
compose  a  macro.  A  macro  is  a  file  of  commands  referenced  by 
the  file  name  in  which  the  statements  are  stored.  For 
example,  if  the  statements  of  Fig  III-l  are  stored  in  file 
COMP,  the  macro  is  executed  by  COMP  DIATOM  DIA  or  COMP.D 
DIATOM  DIA.  In  the  first  example,  the  program  stored  in  the 
tile  DIATOM  (Appendix  F)  is  compiled  and  linked.  The 
executable  module  is  stored  in  the  file  DIA.  In  the  second 
example  DLATOM  is  compiled  using  the  debug  (D)  compiler 
option.  The  source  code  in  the  file  DIATOM  may  either 
contain  all  user  routines,  or  references  to  other  files 
containing  more  source  code.  The  $ADD  facility  causes 
another  file  (filename  =  xxxx)  to  be  inserted  where  $ ADD  xxxx 
occurs  (for  examples,  see  Appendices  B,  D,  F  and  H) .  The 
executable  code  is  run  by  entering  the  name  of  the  file 
containing  it. 

Before  executing  the  program,  the  user  must  attach  all 
input/output  files.  If  the  program  reads  from  logical  unit 
II  and  writes  to  logical  unit  13,  then  they  must  be  attached 
to  the  tries  by  the  AS  (assign)  command  (e.g.  AS  11  =  INFN 
and  AS  13  =  OU'i’FN).  INFN  and  OUTFN  refer  to  the  file  names 
involved.  The  macro  shown  in  Figure  III-2  makes  all 
necessary  assignments  and  executes  the  program  DIATOM  stored 
in  executable  iorm  in  the  file  DIA. 
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III.  Computer  Programs 

This  section  presents  the  computer  code  written  to 
implement  the  theory  outlined  in  Chapter  II.  The  programs 
developed  are  called  DUNHAM,  EMIT,  DIATOM,  and  E' CP ACT.  The 
environment  these  programs  run  in  is  first  discussed.  This 
includes  the  compile,  link,  and  execution  steps.  Then  each 
individual  program  is  discussed,  including  its  purpose, 
cababi 1 i t ios ,  and  major  sub-program  tasks.  Finally,  the 
interrelationships  between  the  programs  is  presented. 
Program  flow  diagrams  and  source  code  are  in  the  Appendicies 
A  -  II.  Appendix  I  is  a  legend  for  the  flow  diagrams. 

Program  Environment 

The  programs  run  on  AFIT's  Harris  800  minicomputer  under 
the  VOS  operating  system.  Each  program  may  run  in  either  a 
batch  or  interactive-terminal  environment.  The  Harris  is  a 
24  bit  machine  with  two  words  used  to  represent  a  real 
variable.  Single  precision  mathematics  is  used. 

The  code  is  written  in  standard  Fortran  77  and  should  be 
portable.  A  few  Harris  utility  routines  are  used,  but  they 
do  not  affect  any  calculation.  Therefore,  they  may  be 
replaced  by  similar  routines  found  on  other  machines  or  left 
out  of  the  programs.  Standard  routines  from  the  IMSL  (11) 
library  are  used.  These  routines  may  also  be  replaced  if  not 
aVd i lable . 
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9-1  CM 


a  wave  function  of  one  electronic  state 

=  a  wave  function  of  the  other  electronic  state 

S  ~  the  matrix  S  of  Fi^.  II-7  and  Eqs  (49)  and 

(62) 


Fi^.  11-11.  The  Franck-Condon  Calculation  in  Matrix  Form 
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The  development  of  Franck-Condon  factors  depends  on  the 
born-Oppenheimer  approximation  (10:199)  that  the  wave 
function  is  separable  into  an  electronic  and  vibrational  wave 
function  4/  =  e  i|r  v  ^  Also,  the  electron  is  assumed  to 
change  states  so  quickly,  compared  to  vibrational  motion, 
that  the  nuclei  have  nearly  the  same  position  and  velocity 
before  and  after  the  transition  (10:199). 

The  inner  product  in  Eq  (75)  is  calculated  in  the  same 
manner  as  the  wave  functions  were  normalized.  In  fact,  the 
same  matricies  can  be  used.  The  difference  arises  since  two 
different  vectors  representing  the  wave  function  are 
involved.  The  problem  illustrated  in  Figure  II-7  is  changed 
to  that  in  Figure  11-11. 


Fig.  11-9.  Franck-Condon  Parabola 
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(49)  illustrated  in  Figure  1 1-7  as  £TS  z.  This  inner  product 
is  constrained  to  be  one  (i.e.,  N  m  ust  be  1)  by  the 
Lagrungian  multiplier  constraint  of  Eq  (63)  as  zrS  z  =  1, 
which  is  used  to  define  the  eigenvalue  problem. 

Cdlculat ion  of  Franck-Condon  Factors 

The  wave  functions  for  the  different  vibrational  states 
of  each  electronic  state  are  now  used  to  calculate  the 
transition  probability  between  two  states  ^ v1  and  ^  v".  This 
probability  is  called  the  Franck-Condon  factor  and  is  defined 
as  (22:119): 

^v'.v"  ~  J  Wdr 

There  are  two  notable  properties  of  the  Franck-Condon 
factor.  First,  the  factors  are  normalized  so  that  (22:120): 


Second,  the  maximum  transition  probabilities  lie  along  a 
parabola— 1  ike  curve  in  the  array  of  factors  (10:196).  This 
is  depicted  in  Figure  11-9.  The  transition  probability  is 
highest  for  the  wave  functions  which  overlap  the  best.  As 
Figure  11-10  illustrates,  there  are  two  wave  functions  Iv" 
for  which  overlap  is  a  maximum  with  +  v  ’ .  This  occurs  for 
all  but  the  v  =  0  wave  function  where  there  is  only  one 
maximum.  This  corresponds  to  the  vertex  of  the  Condon 


parabola . 


Once  recovered  from  y_  and  L,  the  vector  z  contains  the  values 
and  slopes  of  the  function  U|rj  at  each  node  r^.  The  values 
of  Ujrj  are  the  values  of  the  one  dimensional  wave  function. 
These  wave  functions  should  still  be  normalized  it  the 
vectors  were  normalized,  since  the  recovery  processes  were 
unitary.  however,  the  wave  functions  are  checked  for 
nor ma  1  ization  before  using  them  to  calculate  Franck-Condori 
factors.  The  original  three-dimensional  (r,  0,0)  wave 
function  is  not  needed  since  symmetry  is  assumed  and 

is  the  one-dimensional  (r)  wave  function. 

The  eigenvectors  '£  should  be  orthonormal  since  the 
Li  anstor  mations  from  eigenvector  were  unitary  and  is 
orthonormal.  The  normality  of  the  function  U  ^  r  j  ,  as 
described  by  the  vector  z,  is  checked  by  calculating  a 
normal ization  constant  N  from  the  inner  product  of  U(r)  as: 

OO 

,N|2  =  <J0drUn)U(r))  1721 


If  U  ( r )  is  normalized,  N  =  1.  If  U(r)  is  riot  normalized,  it 
can  be  made  so  (U^(r))  by: 


Vrp  "  ""(rj) 


(73) 


The  inner  product  of  U(rj  was  originally  defined  as 


“  £drU(r) 


(74) 


in  the  denominator  of  Eq  (25).  In  matrix  form,  Eq  (74)  is  Eq 


The  matrix  X  is  built  similarly.  Since  Z  =  X  LT,  is 

also 


j 

=  2  X 

k=j-3 


lt 

ik^kj 


(68) 


The  sum  is  again  limited  since  LT  has  only  four  non-zero 
diagonals  also.  Remembering  that  =  Lj^.,  expanding  Eq 
(6b)  and  inverting  yields 


xii=(zl.rxK.i-3)I'.i(.i-3)~xU.i-2)L.i(.i-2)'xi(j-i)I'.i(.i-l))  (69) 

1J  L  •  • 

JJ 

The  eigenvectors  of  Eq  (65)  are  the  energy  levels 
desired,  but  the  eigenvectors  ^  do  not  contain  direct 
information  about  U|r).  This  information  can  be  recovered 
though . 

Since  ^  is  related  to  z  by  ^  =  LT  z,  the  element  y^j  is 


1+3  m 


i+3 


k=i 


kizk  j 


(70) 


Expansion  and  inversion  of  Eq  (70)  yields: 


0 


« 

[  ' 


u 


=  Hz  -  ALLTz 
0  =  L_1H  LT-1LTz  -  ^Ll1LLTz 

i  - 

- 1  T  - 1 

0  =  L  Hi/ 

-1  T-l 

X  =  L  AHLA  (65, 

0  =  X£  -  k}L 


The  decomposition  of  S  is  the  same  used  to  decompose  A  in  the 
previous  section,  and  will  not  be  presented  again.  It  should 
be  noted  though  that  L  is  a  lower  triangular,  banded  matrix. 

The  matrix  X  is  derived  in  the  following  two  step 
process.  First  a  matrix  2  is  found  where  r2_  -  L-1  H.  Then  X 
is  related  to  Z  by  X  =  Z 

The  first  step  is  written  H  =  L  Z  where  the  member 
of  H  is: 


k4:-3ikZkj 


(66) 


The  sum  over  k  is  limited  since  L  has  only  four  non-zero 
diagonals.  Expanding  Eg  (66)  and  inverting  yields 


2i.i~(Hi.i~Li(i-3)Z(i-3)rLi(i-2)2(i-2),rLi(i-l)Z(i- 


L- 

li 


111 


:)  (67) 
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0. 


yields  li  z  -\  S  z  = 

The  eigenvalues  of  Eq  (62)  are  the  system's  expected 
energy  levels.  The  eigenvectors  contain  the  value  of  the 
function  Ujrj  r.rid  its  slope  at  each  node  r^. 

The  solution  of  Eq  (62)  requires  significant  computer 
resources  for  large  matricies  H  and  S.  AFIT  does  not  possess 
computer  code  to  take  full  advantage  of  the  banded  nature  of 
the  matricies  involved.  Over  96%  of  the  members  of  H  and  S 
are  zero,  and  only  half  of  the  4%  non-zero  members  need  be 
involved  in  the  solution  due  to  symmetry. 

A  more  efficient  solution  is  now  presented  that  allows 
the  number  of  operations  to  be  minimized  for  this  problem. 
The  generalized  eigenvalue  problem  H  z  -  \  S  z  =  0  is 
transformed  to  the  standard  eigenvalue  problem  X  ^  -k  %  =  0. 
AFIT  possesses  computer  code  to  efficiently  solve  this 
problem. 

First,  the  Cholesky  decomposition  is  used  on  S  since  it 
is  symmetric,  positive  definite  (28:229).  A  matrix  L  is 
found  such  that  S  =  L  LT.  Then  the  problem  evolves  as 


f ol lows : 


The  extrema  of  E  occur  when  its  first  derivatives  are 


zero,  that  is  when: 


0  =  uE  (59) 

T 

dz 

0  =  Hz  -  ( zTHz)Sz  (60) 

zTSz  (zTSz)2 


0  =  (zTSz)Hz  -  (zTHz)Sz  (61) 

(zTSz)2 

0  =  Hz  -  kSz  (62) 

where  A.  =  (zT  H  z)  /  (zT  S  z)  =  E,  the  energy  eigenvalues.  Eq 
(62)  is  the  generalized  eigenvalue  problem  where  S  is  a 
positive  definite  matrix.  Notice  that  both  H  and  S  are  real, 
band  symmetric  niutricies. 

An  alternative  approach  arrives  at  the  same  generalized 
eigenvalue  problem.  The  method  of  Lagrangian  multipliers  can 
be  used  to  impose  the  normalization  constraint  (z^S  z  -  1  ~ 
0).  Then  the  problem  becomes: 

£  ~  extiz JHz  -  A.(z*Sx  -1)}  (63) 


where  ^  is  a  set  of  multipliers.  Then  taking  the  first 
derivative  equal  to  zero 
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(56) 


(57) 


The  problem  originally  described  by  Eq  (21)  can  now  be 
written  as: 


E  =  exti 


T 

zxH  z 


(5b) 


T 

z  S  z 


where  H  =  T  +  is  the  sum  of  matricies 

describing  the  system's  kinetic  and  potential  energy.  Matrix 
T  is  built  from  the  matricies  representing  I-p  and  the  V  's 
are  built  from  the  matricies  representing  Die  l2a's 
same  manner  as  S  was  built  from  I3  in  Figure  11-7. 
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File  XDIA 


MS 

FR  ALL 
AS  U=&1 
AS  1 2=42 
AS  13=43 
AS  14=IN14 
AS  15=DIARES 
AS  16=DIAP0T 
AS  17=DIAWAV 
AS  99 =$1 3 
AS  o=$Sl3 

PR  Program  DIATOM  is  now  executing  , 
DIA 

PR  ancl  is  now  done. 

ME 

to  use  enters  XDIA  INll  IN12  OUT13 


Fig.  I I 1-2.  Program  Execution  Macro 


Program  DUNHAM 


The  program  DUNHAM  calculates  energy  levels  by  the 
Dunham  Eq  (1)  given  in  Chapter  II.  These  level",  may  be 
needed  to  execute  DIATOM. 


D'JNHAM  accepts  up  to  one  hundred  coefficients  Y^j  where 
i  and  j  range  from  0  to  9.  These  coefficients  are  used  to 
calculate  up  to  676  energy  levels  T(v,J)  (v  =  0  to  26;  J  =  0 
to  26)  using  Eq  (1).  All  levels  may  be  shifted  a  uniform 
amount  so  the  lowest  level  T(0,0)  is  any  value  the  user 
desires.  The  user  controls  the  number  of  energy  levels 
calculated  with  data  elements  VIBLMT  and  ROTLMT.  These  data 


elements  are  the  largest  allowed  values  of  v  and  J.  The 


TiWi  I 


number  of  energy  levels  written  to  the  output  file  is  limited 
by  the  data  element  LVLLMT.  These  levels  are  written  in  the 
format  required  by  DIATOM. 

DUNHAM  is  composed  of  the  sub-program  modules  listed  in 
Table  Ill-l.  The  Harris  routines  (noted  by  a  '*')  are  not 
required  for  energy  calculations. 

The  program  flow  is  presented  in  Appendix  A.  The  source 
code  is  in  Appendix  B.  DUNHAM  first  calls  HDRDUN  which  uses 
BT1ME  and  STIME  to  keep  track  of  CPU  use  and  how  long  the 
program  runs.  Then  HDRDUN  opens  the  output  listing  file 
(unit  13)  and  writes  the  listing  header. 

RDRDUN  is  then  called  to  read  all  input  data  from  unit 
11.  The  input  file  (unit  11)  contains  fixed  form,  key-worded 
records  (card  images)  as  shown  in  Figure  I1I-3.  Data  records 
are  identified  by  a  ’>'  in  column  1.  All  other  records  are 
ignored  as  user  comments.  Each  data  record  must  contain  a 
valid  key  word  in  columns  2  to  4.  The  valid  key  words  and 
related  data  elements  are  given  in  Table  III-2.  The  value  of 
all  data  elements  except  HDR  must  be  in  columns  6  to  20.  The 
character  string  used  for  HDR  must  be  in  columns  6-35.  All 
real  data  elements  must  contain  a  decimal  point.  For 
example,  ’>Y10=100'  is  invalid,  whereas  ’>Y10=100.0'  is 
valid.  Integer  data  elements  may  not  have  a  decimal. 
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Table  III-l 


Program  DUNHAM  Modules 


MODULE 

TYPE 

DESCRIPTION 

DUNHAM 

Main 

The  main  module  -  calculates  the 
energy  levels. 

HDRDUN 

Subroutine 

Opens  output  listing  file. 

BT11VIE  * 

it 

Starts  CPU  use  statistics. 

STIME  * 

If 

Starts  wall  clock  use  statistics 

DATE  * 

it 

Returns  the  current  date. 

TIME  * 

ll 

Returns  the  current  time. 

USERNO* 

II 

Returns  the  user's  ID. 

RDRDUN 

II 

Controls  all  input. 

TRLDON 

ll 

Closes  the  output  listing  file. 

ETIME  * 

II 

Stops  CPU  use  statistics. 

WTIME  * 

ll 

Stops  wall  clock  use  statistics. 

*  Harris  Routines  (9) 


RliO  X 


15  OCT  64 


w 


■»***  JEFF  +  CONRAD’S  DATA  -- 

>R0T=0 

>VXB=25 

>LVL=25 

>HDR=PBO  X  state 
>DEQ-0. O 
>Y00“0 

>Y10«722. 687 
>Y20«  -3. 613 


Data  Elements 


Key  Words 


Fig.  III-3 .  Program  DUNHAM  Input  File 


Table  III-2 


Program  DUNHAM  Key  Words  for  Input  Data 


Key  Word 

Data  Element 
Type 

Data  Element 

VIB 

Integer 

Vibrational  quantum  number  -  upper 
limit . 

HOT 

Integer 

Rotational  quantum  number  -  upper 
limit. 

LVL 

Integer 

Total  Number  of  energy  levels  to  be 
written  to  the  output  file. 

HDR 

Character 

A  label/comment  written  to  the 
output  file. 

DEQ 

Real 

The  energy  value  of  the  dissocia¬ 
tion  limit. 

Yij 

Real 

A  Dunham  coefficient  (e.g.  Y10  is 
the  Key  word  for  Y^0)  • 
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After  returning  from  RDRDUN/  DUNHAM  calculates  each  term 
of  Eq  (1)  and  adds  it  to  the  variable  SUM.  Once  all  values 
of  Eq  (1)  are  collected,  the  value  of  the  energy  level  is 
transferred  to  the  variable  T.  After  all  allowed  energy 
levels  T  have  been  calculated,  they  are  shifted  so  the 
oissociation  energy  level  value  is  the  value  input  on  the 
'DEQ'  record.  The  shifted  energy  levels  are  then  written  to 
the  output  file  (unit  14)  in  the  format  required  for  DIATOM. 

Finally,  TRLDUN  prints  the  run  statistics  (CPU  use  and 
execution  time)  and  closes  the  output  listing. 

Program  EFIT 

The  program  EFIT  uses  measured  transition  lines  in  a 
weighted  least  squares  calculation  of  observed  energy  levels. 
These  energy  levels  are  used  as  input  to  DIATOM. 

EFIT  accepts  62,500  transition  lines  involving  250 
energy  levels.  A  maximum  of  ten  electronic  states  may  be 
entered  with  up  to  25  levels  for  each  state.  These  levels  do 
not  have  to  be  the  first  25  (v  =  0  to  24),  but  can  be  any 
level  as  long  as  they  are  properly  labelled  on  input.  The 
user  can  shift  all  levels  so  that  the  lowest  level  of  the 
lowest  state  has  a  specified  value.  Each  level  calculated  is 
written  to  an  output  file  in  the  format  required  by  DIATOM. 

EFIT  is  composed  of  thtr  sub-program  modules  presented  in 


Table  III-3.  The  program  flow  is  depicted  in  Appendix  C  and 
the  source  code  is  in  Appendix  D. 


Table  II1-3 


Program  EFIT  Modules 

MO  DULE 

TYPE 

DESCRIPTION 

EFIT 

Main 

The  main  module  -  does  the  least 
squares  fit. 

HDREFT 

Subroutine 

Opens  the  output  listing  file. 

BTIIVIE  * 

II 

Starts  CPU  use  statistics. 

STIME  * 

H 

Starts  wall  clock  use  statistics. 

DATE  * 

It 

Returns  the  current  date. 

TIME  * 

II 

Returns  the  current  time. 

USERNO* 

II 

Returns  the  user's  ID. 

RDREFT 

II 

Controls  all  input. 

LVLEFT 

If 

Calculates  the  absolute  energy- 
level  number  for  an  energy  level 
relative  to  an  electronic  state. 

TRLEFT 

II 

Closes  the  output  listing  file. 

ETIME  * 

II 

Stops  CPU  use  statistics. 

WTIME  * 

II 

Stops  wall  clock  use  statistics 

*  Harris  Routines  (9) 


EFIT  first  calls  HDREFT  which  performs  the  same 
functions  described  for  HDRDUN  in  program  DUNHAM. 

Then  RDREFT  opens  and  reads  all  records  in  the  input 
file  (unit  11).  The  input  technique  is  the  same  described 
for  RDRDUN  of  DUNHAM.  The  key  words  are  six  characters  long 
and  are  in  colums  2-7  of  each  record.  Table  I l 1-4  gives  the 
acceptable  key  words  and  related  data  elements. 

Figure  1 1 1  —  4  is  an  example  of  a  valid  input  file.  In 
this  example,  the  first  state  (STAT01)  is  the  X  or  ground 
electronic  state,  the  second  (STAT02)  the  a  or  first 
electronically  excited  state,  and  so  on.  Notice  that  data 
for  the  X  state  vibrational  levels  v  =  0  to  8  involved  in  the 
following  records  are  indicated  by  the  record  keyworded 
LVLS01.  Also  notice  that  the  b  state  has  only  one  level 
involved  (v  =  4).  Missing  levels  are  allowed.  The  program 
does  not  require  every  level  to  be  represented  by  the  data  in 
the  range  v  =  (lowest  value)  to  v  =  (highest  value)  for  each 
state.  The  record  keyworded  SHIFTS  contains  the  value  of  the 
lowest  level  of  the  lowest  state  on  output.  All  other  levels 
will  be  shifted  so  this  occurs.  Transition  data  is  given  by 
the  remaining  records  key-worded  ABBCDD.  These  records 
contain  the  transition  line  observed  (cm-1)  and  optionally 
the  least  squares  weighting  factor  (0  w  _<  1.0)  separated  by 
a  for  the  transition  from  vibrational  level  BB  of 
electronic  state  A  to  level  DD  of  state  C.  If  the  weighting 
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Table  III-4 


Program  EFIT  Key  Words  for  Input  Data 
DATA  ELEMENT 


KEY  WORD 

TYPE 

DATA  ELEMENT 

STATxx 

Character 

A  single  character  symbol  for 
the  electronic  state  xx;  e.g., 
STAT01 =x  assigns  the  symbol 
"x"  to  the  first  state. 

LVLSxx 

Integer 

A  series  of  integer  numbers 
(1  or  2  digit)  indicating  which 
energy  levels  V=?  are  involved 
in  the  least  squares  fit  for 
electronic  state  xx  symbolized 
by  STATxx. 

SHIFTS 

Real 

The  value  of  the  lowest  level 
of  the  lowest  state. 

ABBCDD  ljw 

Real 

Two  data  elements  -  the  value 
(1)  of  the  transition  observed 

between  vibrational  level  BB 
of  electronic  state  A  and  level 
DD  of  state  C  and  its  weighting 
factor  (W^O  to  1.0).  The  two 
values  are  separated  by  a  "t". 


51 


Fig.  III-4.  Program  EFIT  Input  File 


factor  is  missing  on  input  it  defaults  to  1.0,  the  highest 
weighting  allowed.  The  order  of  the  keywords  is  important 


and  should  match  that  given  by  Figure  III-4.  For  example, 
the  LVLS01  record  can  not  occur  before  STAT01  record. 

Control  now  passes  back  to  the  main  module.  EFIT  builds 
the  matrix  A  (Figure  1 1-3)  of  the  least  squares  equation  A  x 
-  b.  A  is  stored  in  full  matrix  storage  mode  such  that  A(1,J) 
corresponds  to  Axj.  The  matrix  b  is  built  next.  Then  the 
problem  is  changed  to  =  b  where  A  =  L  L1'  and  y  =  LT  x  by 
the  Cholesky  decomposition  Eq  (10).  The  solution,  is 
found  through  Eq  (12)  and  transformed  to  x  by  Eq  (15)  after  a 
consistency  check  is  made  (Eq  (13)).  The  matrix  x  contains 
the  unshifted  energy  levels.  Then  all  levels  are  shifted  so 
the  lowest  level  is  at  0.  The  levels  are  then  written  to  the 
output  listing.  Finally,  they  are  shifted  again  so  that  the 
lowest  level  is  at  the  value  input  by  the  SHIFTS  records. 
These  final  energy  level  values  are  written  to  an  output  file 
suitable  for  DIATOM. 

A  possible  source  of  confusion  is  the  energy  level 
labelling  used  by  EFIT.  EFIT  gives  each  level  an  absolute 
level  number,  while  also  maintaining  a  level  number  relative 
to  the  electronic  state  for  the  user's  convenience.  This  is 
illustrated  in  Table  II 1-5.  EFIT  uses  absolute  level  numbers 
ifi  array  references  and  calculations. 
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Table  II1-5 


Relative  and  Absolute  Energy  Level  Numbering 

STATE 
NUMBER 

01 

02 

03 
04 

Program  D 

The  program  DIATOM  finds  a  set  of  parameters  for  the 
potential  energy  model  that  best  fits  the  observed  energy 
levels  in  a  weighted  least  squares  sense.  Then  DIATOM  finds 
the  normalized  wuvt'  functions  for  the  calculated  energy 
levels.  These  wave  functions  are  needed  for  the  program 
FCF ACT . 

The  execution  of  DIATOM  may  be  modified  by  keyworded 
records  'RM'  and  'KW'.  If  RM  -  1,  DIATOM  finds  the  potential 
energy  parameters  that  best  fit  the  energy  levels.  If  RW  = 
1,  DIATOM  calculates  the  wave  functions. 

DIATOM  accepts  up  to  25  observed  vibrational  energy 
levels  for  a  specific  electronic  state.  However,  these 
levels  must  be  consecutive  and  start  at  v  =  0.  Missing 


STATE 

SYMBOL 


b 

A 


RELATIVE 

NUMBER 

1 

2 

3 

1 

2 


0 

1 

2 

3 

4 
8 


ABSOLUTE 

NUMBER 

1 

2 

3 

4 

5 


? 

8 

9 

10 

11 

12 


I  ATOM 
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levels  may  be  marked  by  assigning  it  a  value  of  zero.  Then 
no  comparison  is  made  by  DIATOM  to  the  calculated  value.  The 
grid  used  by  DIATOM  may  have  up  to  100  steps  (101  nodes). 
The  potential  energy  model  may  use  up  to  ten  constants  and 
ten  parameters.  Each  parameter  is  associated  with  an  upper 
and  lower  limit. 

DIATOM  is  composed  of  the  modules  listed  in  Table  III-6. 
The  program  flow  is  presented  in  Appendix  E  and  the  source 
code  is  in  Appendix  F. 

DIATOM  first  calls  HEADER  which  performs  the  same 
functions  HDRDUN  of  program  DUNHAM  performed. 

Then  READER  opens  and  reads  all  records  from  the  input 
files  (unit  11,  12,  and  14).  The  input  technique  described 
for  RDRDUN  of  program  DUNHAM  applies.  Unit  11  contains  data 
that  control  the  execution  of  DIATOM.  Valid  key  words  and 
data  elements  are  listed  in  Table  III-7.  The  keywords  must 
be  in  columns  2-3.  An  example  of  a  valid  unit  11  input  file 
is  given  in  Figure  II1-5.  Records  key  worded  'Cx'  contain 
constants  used  by  the  potential  energy  model.  Records  key 
worded  'Lx',  'Px',  and  'Ux'  contain,  respectively,  values  of 
the  lower  limit,  initial  value,  and  upper  limit  of  a 
potential  energy  parameter. 


Unit  12  contains  data  related  to  the  observed  energy 
levels.  Valid  key  words  and  data  elements  are  given  in  Table 
Ill-b,  and  an  example  of  a  valid  unit  12  input  file  is 


Program  DIATOM  Module:; 


MO DULE 

TYPE 

diatom 

Main 

HEADER 

Subroutine 

HTIME 

* 

»< 

STIME 

* 

II 

DATE 

* 

H 

TIME 

* 

II 

USERNO 

* 

II 

READER 

II 

MiNUM 

II 

FUN 

Pune  tion 

POTENT 

Subrou  tine 

EICEN 

II 

CKTL 

ll 

FOLDZ 

ii 

FOLDX 

ll 

EHOUSS 

tt 

ll 

EQRT1S 

it 

it 

E'QKTzS 

it 

ll 

EllO  HRS 

tt 

M 

UNFOLD 

ll 

PUTRND 

ll 

WaVE 

li 

ICSCCU 

it 

ll 

No KMAL 

H 

VMULQP 

it 

ii 

VI.lULFF 

tt 

ii 

PUT WAV 

ll 

0  UTl’UT 

M 

CETPUT 

II 

PUNTER 

II 

PLTRES 

II 

PLTPOT 

II 

tkailr 

II 

ET I  ME 

Mt 

II 

WTIME 

* 

II 

*  Harris  Hou tinea  (9) 
ft  IMSL  Routines  (11) 


_ DESCRIPTION _ 

The  main  program  module. 

Opens  the  output  listing  file. 
Starts  CPU  use  statistics. 

Starts  wall  clock  use  statistics. 
Returns  the  current  date. 

Returns  the  current  time . 

Returns  the  user’s  ID. 

Controls  all  input. 

Finds  the  minimum  of  FUN. 

The  sum  of  weighted  residuals 
squared. 

Returns  potential  energy  value 
and  slope. 

Solves  H_v  -  S v  :::  0 . 

Cholesky  decomposition  S  =  LL1 . 
Finds  Z  where  H  -  LZ . 

Finds  X  where  X  -  Zh'"'^ . 

Changes  X  to  tridiagonal  foi’in. 
Finds  eigenvalues  of  tridi¬ 
agonal  X. 

Finds  eigenvalues  and  eigen¬ 
vectors  of  tridiagonal  X. 

Finds  eigenvectors  of  X  from 
those  of  tridiagonal  X. 

Finds  eigenvectors  v  from  those 
of  X. 

Stores  random  number  genera¬ 
tor  values. 

Finds  the  wave  functions. 

Fits  a  cubic  spline  to  a  function 
Normalizes  the  wave  functions. 
Matrix  multiplication. 

Matrix  multiplication. 

Stores  wave  functions  in  a 
file . 

Writes  output  listing  and 
da  ta . 

Calculates  1000  potential  ener¬ 
gy  values. 

Creates  the  listing. 

Creates  residual  plot  file. 
Creates  potential  plot  file. 
Closes  the  output  listing  file. 
Stops  CPU  use  statistics. 

Stops  wall  clock  use  statis¬ 
tics. 


Table  III-7 


Program  DIATOM  Key  Words  for  Input  Data  -  Unit  11 


KEY  WORD 

TYPE 

DATA  ELEMENT 

BG 

Real 

The  leftmost  (beginning)  mode  of 
the  grid. 

EN 

H 

The  right  most  (ending)  inode  of 
the  grid. 

NE 

In xeger 

The  number  of  grid  elements 
( steps) . 

IS 

It 

The  number  of  steps  MINUM  may 
take . 

IP 

ll 

The  modulus  of  the  steps  at 
which  MINUM  printing  is  required 
If  negative  no  information  is 
printed. 

JR 

H 

Random  step  frequency  (MINUM) . 

JG 

ll 

Gradient  step  frequency  (MINUM). 

JA 

M 

Average  step  frequency  (MINUM). 

JJ 

H 

Jump  step  frequency  (MINUM). 

FR 

il 

A  flag  for  creating  residual 
plot  file  (FR-'l  yes;  FR=-0  no). 

FP 

•  I 

A  flag  for  creating  potential. 

klvl 

Integer 

A  flag  controlling  execution  of 
MINUM.  If  RM-'l,  MINUM  is  used. 
If  RM~0 ,  MINUM  is  not  used. 

RW 

tl 

A  flag  controlling  execution  of 
WAVE.  If  RW~1,  WAVE  is  used. 

If  RW--0,  WAVE  is  not  used. 

It  B 

Real 

The  value  of  Planck's  Constant  h 

MU 

il 

The  molecule’s  reduced  mass. 

Gx 

li 

"x"  =  0  to  9«  The  value  of  one 
of  10  constants  available  to  the 
potential  energy  model. 

Px 

ii 

"x"  :  0  to  9*  The  initial  value 
of  one  of  ten  paramenters  avail¬ 
able  to  the  potential  model. 
MINUM  changes  these. 

Lx 

il 

Lower  limit  of  Px. 

Ux 

li 

Upper  limit  of  Px. 
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CLCFCF  builds  the  matricies  of  Figure  11-11  and  multiplies 
them  in  a  manner  similar  to  that  used  by  NORMAL  of  DIATOM. 
The  matrix  product  is  the  inner  product  between  the 
normalized  wave  functions.  This  product  is  squared  and 
stored  in  the  two  dimensional  array  FACTOR  as  the  Franck- 
Condon  factor. 

OUTFCF  prints  the  Franck-Condon  table  in  two  sections 
after  the  first  loop  is  done.  Then,  TRLFCF  stops  the  run 
statistics  and  closes  the  output  listing  file  (unit  13) . 

Program  Relationships 

The  programs  DUNHAM,  EFIT,  DIATOM  and  FCFACT  are  related 
as  shown  in  Figure  III-ll.  The  user  may  use  DUNHAM,  EFIT,  or 
some  other  source  of  energy  level  data  to  build  the  input 
file  (unit  12)  for  DIATOM.  DIATOM  in  turn  builds  one  of  two 
input  files  for  FCFACT.  The  second  FCFACT  input  file  is 


built  by  running  DIATOM  for  a  different  electronic  state  (a 
new  set  of  energy  levels  and  potential  parameters) . 


CO  00 


>LABL«Singl«  Harmonic 
>S  1  0.000000 

>S  2  0.000000 

>S  3  0.000000 

>S  4  -O.  3638484E -02 
; 3  5  0.000000 

>  S  6  -O. 363S484E -02 

>S  7  0.000000 

>  S  8  -O. 3638484E  02 

•S  9  0.000000 

>S  lO  -O. 3638484E -02 
'5  11  0.000000 

>S  12  -O. 3638484E -02 
S  13  0.000000 

>5  14  -0.3638484E  02 
S  It.  0.000000 
' S  16  -0. 3o38484E  02 
S  17  0.000000 

>S  18  -O. 3638484E  -02 
>S  19  0.000000 

/ 3  20  -O. 3638484E- 02 
-S  21  0.000000 

-S  22  -  0. 3o384tl4E- 02 
S  23  0.000000 


Gs»c l 1 1  fetor 
0. 000000 
0.000000 
O. 4406163E-01 
-0. 267S800E-03 
0.4408163E-01 
-0. 2876600E-03 
O. 44081&3E-01 
-O. 2878600E-03 
O. 44081o3E-01 
-O. 2878800E-03 
O. 4408163E-01 
-O. 2878800E-03 
O. 44G8163E-01 
-O. 2878800E -03 
O. 4408163E-01 
-0. 2378800E-03 
O. 4408163E-01 
-0. 2878800E-03 
O. 4408163E-01 
-0. 2378800E-03 
0. 44O8163E-01 
-0. 2378800E-03 
O. 4408163E-01 


'8  24  -0.3638484E  02  -0 . 237S300E -03 
S  2b  0.000000  O. 4408163E-01 

•S  26  -0.3638484E  02  -O . 2878a00E -03 


27  0.000000  O. 4408163E-01 

2d  -0. 3o38484E- 02  -O . 2878800E -03 


0. 000000 
O. 6157434E-02 
0. 3638484E-02 
0.000000 
O. 3638484E-02 
O. OOOOOO 
O. 3638484E-02 
O. OOOOOO 
0. 3638484E-02 
0.000000 
0. 3638484E-02 
0.000000 
O. 3638484E-02 
O. OOOOOO 
0. 3638484E-02 
O. OOOOOO 
0. 3638484E-02 
0.000000 
O. 36  38484E - 02 
O. OOOOOO 
0. 3638484E-02 
O . OOOOOO 
0. 3cj38484E-02 
0.000000 
0. 3o38484E-02 
0. OOOOOO 
O. 3638484E-02 
O. OOOOOO 


Matrix  S  Data 


■  3  69 
70 
S  71 
>S  72 
O 

>  O 

o 

•>  o 

>  o 

•)  o 

o 


0.000000  o. 

-0. 3638484E-02  -O. 

0.000000  0. 

-0. 3638484E-02  -O. 


1  -O.  1293752E- 10 

2  -O. 2054230E-09 

3  -0. 846 1519E- 10 

4  -O. 1223209E-08 
3  -0. 1391 142E-09 

6  -O. 3602338E-03 

7  -0. 4&73904E -09 


4408163E-01 
2878800E-03 
44  08 1 63E -  0 1 
2878800E-03 


0. 3638484E-02 
O. OOOOOO 
O. 3638404E-02 
-O. 6 1 374  34E - 02 


Eigenvector  Data 


O. 1273469 

O. 3838401E-03 

0.2346939 

O. 767680 IE -03 

0.2546939 

O. 7676801E-03 

O. 2546939 

O. 767660 1 E - 03 

0.2546939 

0.  /'676801E-G5 

0. 254o939 

O. 76768018-03 

0.2546939 

0. 767o801E-03 

0.2546939 

0. 7676801E-03 

0.2546939 

O. 7676801E-03 

0.2546939 

O. 767680 IE- 03 

0.2546939 

O. 7676801E-03 

0.2546939 

O. 7o7o80 1 E-03 

0.2546939 

0. 7676801E-03 

0.2346939 

0. 76768018-03 


0.2546939 
O. 7676801E-03 
0. 1273469 
O. 3838401 E - 03 


Fig.  III-10.  Program  FCFACT  Input  File 
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Valid  key  words  (columns  2-5)  and  data  elements  are 
listed  in  Table  III-ll.  Ari  example  of  a  valid  input  file  for 
units  11  and  12  is  given  in  Figure  1 11-10.  The  'LABL'  record 
contains  a  72  character  label  of  the  wave  function  set.  All 
records  with  a  "S"  in  column  2  contain  data  belonging  to  the 
matrix  S  of  zj  S  zj_.  The  "SXXX"  record  contains  the  "XXX" 
row  of  matrix  S  in  band  symmetric  form.  The  key  word  of  the 
remaining  "XXXX"  records  is  the  energy  level  number.  Columns 
6-9  contain  a  sequential  record  number  relative  to  the  energy 
level.  Columns  10-24  contain  the  value  of  the  eigenvector 
corresponding  to  the  number  in  columns  6-9. 

Next,  FCFACT  sets  LUNIT  =  12  and  calls  RDRFCF  to  read 
all  data  concerning  the  upper  (v')  state. 

The  Franck-Condon  factors  are  calculated  by  CLCFCF  which 
is  called  for  each  combination  of  wave  functions. 

Table  III-ll 

Program  FCFACT  Key  Words  for  Input  Data 


KEY  WORD 

TYPE 

DATA  ELEMENT 

LABL 

Character 

A  72  character  label  of  the 
wave  function  set. 

"Sxxx" 

Real 

4  values  belonging  to  the 
"xxx"  row  of  matrix  S 
in  band  symmetric  form. 

(e.g.  "SOU"  contains 

S(  11 , 1) »  S(ll,2),  S(ll, 3) , 
and  S(ll,4) ) . 

"XXXX" 

Real 

"XXXX"  =  "1"  to  "25".  The 
data  for  energy  level  "XXXX 
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Table  111-10 


MODULE 

Program  FCFACT  Modules 

TYPE 

FCFACT 

Main 

The  main  module. 

HDRFCF 

Subroutine 

Opens  the  output  listing  file. 

BTIME 

* 

It 

Starts  CPU  use  statistics. 

STIME 

* 

It 

Starts  wall  time  use  statis¬ 
tics  . 

DATE 

* 

II 

Returns  current  date. 

USERNO 

* 

II 

Returns  the  user's  ID. 

RDRFCF 

It 

Reads  one  set  of  wave 
functions . 

CLCFCF 

tt 

Calculates  the  FCF  for  two 
wave  functions. 

OUTFGF 

It 

Prints  the  FCF  table. 

TRLFCF 

It 

Closes  the  output  listing 
file. 

ETIME 

* 

tt 

Stops  CPU  use  statistics. 

WTIME 

* 

M 

Stops  wall  time  use  statistics 

* 

Harris  Routines 

(9) 

67 


Finally,  TRAILR  stops  the  run  statistics  and  closes  the 


output  file. 

Program  FCFACT 

Program  FCFACT  computes  the  inner  product  given  by  Eq 
(75),  (Franck-Condon  factors)  for  two  sets  of  wave  functions. 
Each  wave  function  set  is  the  output  of  one  run  of  DIATOM. 

FCFACT  expects  both  input  files  (units  11  and  12)  to 
have  identical  matricies  S  and  25  (v  =  0  to  v  =  24)  sets  of 
wave  function  data.  Both  sets  must  be  related  to  an 
identical  grid  in  DIATOM.  FCFACT  builds  a  25  by  25  table  of 
Franck-Condon  factors. 

FCFACT  is  composed  of  the  sub-program  modules  listed  in 
Table  III-10.  The  program  flow  is  presented  in  Appendix  G, 
and  the  source  code  is  in  Appendix  H. 

FCFACT  first  calls  HDRFCF  which  performs  as  HDRDUN  of 
DUNHAM  did. 

Then  FCFACT  opens  units  11  and  12,  sets  the  variable 
LUNIT  =  11,  and  calls  RDRFCF.  This  causes  RDRFCF  to  read  all 
data  concerning  the  lower  ( v" )  state. 


where  wi  is  the  weighting  factor  (0  <  w  <  1)  associated  with 


Control  is  passed  to  PUTRND  after  MINUM  finishes. 
PUT  RN  D  stores  the  current  values  of  the  random  number 
generator  seeds  for  next  time  DIATOM  is  executed. 

Next,  WAVE  finds  the  normalized  wave  functions 
(eigenvectors)  for  each  calculated  energy  level  (eigenvalue). 
First,  the  main  module  sets  JOBN  =  1  before  WAVE  is  called. 
Then  WAVE  calls  FUN  using  the  "best"  parameters  found  by 
MINUM.  This  causes  the  eigenvectors  of  the  standard 
eigenvalue  problem  (Eq  (65))  to  be  computed.  WAVE  also 
writes  the  matrix  S  to  unit  17  for  the  program  FCFACT. 

Then  NORMAL  uses  IMSL  routines  VMULQF  and  VMULFF  are 
used  to  multiply  the  matricies  of  Eq  (77),  and  the 
normalization  factor  N  (Eq  (81))  is  computed.  Then  the 
normalized  value  of  the  wave  function  is  computed.  PUTWAV 
then  stores  the  eigenvectors  (wave  functions)  in  unit  17 
using  the  format  required  by  program  FCFACT. 

Then,  OUTPUT  calls  GETPOT  to  calculate  the  value  of  the 
potential  energy  model  at  1000  points  along  the  grid.  These 
values  can  be  used  for  plotting  the  potential.  PRNTER  writes 
the  potential  values,  parameters,  and  energy  data  to  the 
output  listing.  OUTPUT  also  uses  PLTRES  and  PLTPOT  if  DIATOM 
was  directed  to  create  plot  files  of  the  residual  (unit  15) 


and  potential  (unit  16)  data  (FR  =  1;  FP  =  1). 


Then  the  generalized  eigenvalue  problem  Eq  (62)  is 


solved  when  FUN  calls  EIGEN.  EIGEN  first  calls  GETL  to 
decompose  the  matrix  S  by  Eq  (10).  Then  EIGEN  uses  FOLDZ  and 
FOLDX  to  implement  Eq's  (67)  and  (69).  This  is  how  the 
generalized  eigenvalue  problem  is  transformed  into  the 
standard  eigenvalue  problem  Eq  (65).  The  symmetric  matrix  X 
of  the  standard  eigenvalue  problem  is  changed  to  a  symmetric 
tridiagonal  matrix  T  by  IMSL  routine  EHOUSS  (11)  using 
Householder's  transformation  (13).  Then  the  lowest  25 
eigenvectors  of  T  are  found  by  IMSL  routine  EQRT1S  (11)  using 
a  QR  transformation  with  a  Newton  shift  (18).  If  JOBN  =  1, 
the  eigenvectors  of  T  are  computed  by  the  QR  method  (1)  using 
IMSL  routine  EQRT2S  (11).  The  eigenvectors  of  X  are  found 
from  those  of  T  by  IMSL  routine  EHOBKS  (11),  (13). 

Next,  UNFOLD  recovers  the  eigenvectors  of  the 
generalized  eigenvalue  problem  Eq  (62)  from  those  found  by 
FUN. 

FUN  then  compares  the  eigenvalues  v  and  the  observed 
energy  levels  Ev  by  computing  a  residual  (variable  RESID) : 


RESID(I)  = 


( X. I  -  E-j.)  for  E-j-  0 


(78) 


0 


for  Ej  =  0 


Then  the  weighted  sum  the  residuals  is  computed  as  FUN: 


2  5 

FUN  =  iZ  (w?  x  RESID(i)2) 
i=l  1 
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(79) 


Matrix  H 
of  order 

3 


11 

hl2 

h13 

hll 

21 

h22 

h23 

h21 

31 

h32 

h33 

h22 

mm 
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32 

33 


Full 
Storage 
Mo  de 


Symme  ti'ie 
Storage 
Mo  de 


Fig.  III-8.  Symmetric  Matrix  Storage  Mode 


Matrix  S 
of  order 


5 


311 

y12 

0 

0 

0 

0 

sll 

S21 

c? 

22 

s23 

0 

0 

S21 

S22 

0 

c* 

32 

s33 

s34 

0 

- 

s32 

S33 

0 

0 

s43 

s44 

y45 

s43 

y44 

0 

0 

0 

s34 

°55 

s5't 

s55 

Full 
S to  rage 
Mode 


Band 
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Storage 
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Fig.  III-y.  Band  Symmetric  Matrix  Storage  Mode 
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The  variable  JOHN  controls  calculation  of  eigenvalues 
(JOBN  =  0)  or  eigenvectors  (JOBN  =  1).  DIATOM  initially  sets 
JOBN  =  0  since  the  eigenvectors  are  not  needed  until  the  best 
parameter  set  has  been  found  by  MINUM.  MINUM  is  a  routine 
written  by  Pearson  and  Williams  (15)  that  computes  the 
minimum  of  a  real  function  FUN.  FUN  is  a  function  of  ten 
potential  energy  parameters  which  must  be  scaled  so  each 
range  from  0  to  1.0.  MINUM  may  be  replaced  by  any  non-linear 
minimization  routine.  MINUM  takes  random  and  statistically 
derived  "best  guess"  steps  and  jumps  through  FUN's  parameter 
space . 

The  function  FUN  first  builds  the  matricies  H  and  S  of 
Eq  (58).  FUN  calls  POTENT  once  for  each  grid  element  in  this 
process.  The  matrix  H  is  stored  in  symmetric  mode  as  a  one 
dimensional  array  of  length  n(n  +  l)/2  where  n  is  the  order 
of  the  matrix.  The  matrix  element  is  the  k*-*1  element  of 
the  array  where  k  =  (i(i  -  l)/2)  +  j.  Due  to  symmetry,  only 
the  lower  triangle  of  H  (i  j)  is  stored.  This  is 
illustrated  in  Figure  1II-8.  Matrix  S  is  stored  in  band 
symmetric  mode,  illustrated  in  Figure  III-9,  in  a  two- 
dimensional  array.  Only  the  elements  on  the  main  and  sub¬ 
diagonals  are  stored.  The  matrix  element  S^j  (i  =  1  to  n?  j 
=  i,  (i  -  1),  (i  -  2),  (i  -  3))  is  stored  in  the  k,  1^ 
element  of  the  array  where  k  =  i;  1  =  4  -  (i  -  j);  and  j  =  (k 


Table  II1-9 


Program  DIATOM  Key  Word 
Input  Data  -  Unit  14 


KEY  WORD 

TYPE 

DATA  ELEMENT 

IU 

Integer  *6 

Starting  random 
number  seed. 

IX 

II 

Fixed  random 
number  multi¬ 
plier  (ix=131075) 

*#*##<HW#*#*#***###*##*******###*###*#*##***#**#*###***#*********#*## 
*•*#  Seeds  far  a  random  number  generator  in  NINUM 

**#  called  by  program  DIATOM 

frfl**************^-***********!!****************************#******#**#* 
>IU»  5 1 31 34068332b 7 
>IX  =  131075 


Fig.  III-?.  Program  DIATOM  Input  File  -  Unit  14 


*  Single  Harmon ic  Oscillator 

*  24  Aug  84 


>LBL-S ingle 
>V00=1 . 0 
>  V01 =3 . 0 
>V02=5. 0 
> V03-0 . 0" 
>V04=9. O 
>V03=1  1 . 0  i 
>1/06=1  3 .  o  ; 
>vo7=i;>.o  ; 
>V08=  17.0  i 
>009=19.0  I 
>010=21.0  I 


Harmonic  Oscillator 


The  v=3  level  is  not  observed 
and  must  be  entered  as  0. 


The  v=8  level  is  weighted 
6  tenths  of  the  other 
Itvels  (v=0  to  4)  which 
default  to  1.0. 


Fig.  III-6.  Program  DIATOM  Input  File  -  Unit  12 


are  entered  starting  a  v  =  0  to  4  and  v  =  6  to  13,  then  14 
records  'Vxx'  must  be  entered  'xx‘  =  *00’  to  ’13'.  Record 
1  VO 5 '  must  be  zero  since  no  level  was  observed  for  v  =  5. 
This  level  will  not  be  used  in  the  energy  level  fits  and 
serves  only  as  a  place  holder. 

Unit  14  contains  values  used  as  seed  numbers  for  MINUM's 
random  number  generator.  DIATOM  replaces  the  seed  values 
with  new  values  every  run.  Therefore,  the  user  need  only 
supply  this  data  when  the  program  is  first  put  on  the 
computer.  At  this  time  both  may  be  set  to  131075.  Valid  key 
words  must  be  in  columns  2-3  and  are  listed  in  Table  III-9. 


An  example  of  an  input  file  is  given  in  Figure  III-7. 


Table  111-8 


Program  DIATOM  Key  Words  for  Input  Data  -  Unit  12 


KEY  WO  RD 

TYPE 

DATA  ELEMENT 

LBL 

Character 

A  72  character  label 
used  for  the  wave 
function  output  file. 

Vxx  e;w 

Real 

"xx"  =  "01"  to  "10". 
Two  data  elements. 

First  is  the  value 
of  the  observed  energy 
(e)  level  xx.  Se¬ 
cond  is  a  weight¬ 
ing  factor  ( w  —  0  to 
1.0)  applied  to  that 
level  separated  by 
a  "  ;  "  . 


presented  in  Figure  III-6.  Keywords  must  be  in  columns  2-4. 
The  character  string  contained  in  the  'LBL'  record  is  used  to 
label  the  wave  function  output  file  (unit  17).  bach  'Vxx' 
records  contains  the  observed  energy  value  and,  optionally,  a 
weighting  factor  (0  _<  w  <_  1.0)  for  the  energy  level  v  -  'xx'. 
The  units  of  the  energy  level  must  be  the  same  used  for  the 
values  of  Planck's  constant  fi  (record  'lib1  in  unit  11)  and  the 
system's  reduced  mass  H  (record  'MU'  in  unit  11).  The  energy 
level  value  and  the  weighting  factor  are  separated  by  a  1 ; 1 . 
The  weighting  factor  defaults  to  1.0  if  it  is  not  entered. 
The  'Vxx'  records  must  be  consecutive  in  'xx'.  If  13  levels 
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MINUM  is  not  allowed  to  run. 
WAVE  is  allowed  to  run. 

Potential  Energy  parameter  1 
starts  at  2.0  and  may  vary 
from  1.99  to  2.01 


Fig.  III-5.  Program  DIATOM  Input  File  -  Unit  11 


IV.  Hesults  and  Discussion 
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This  section  consists  of  three  parts.  The  first  two 
parts  present  the  validation  of  the  programs  DIATOM  and 
FCFACT.  Then  a  CPU  use  benchmark  is  presented  of  the  program 
DIATOM  in  the  third  part.  The  programs  DUNHAM  and  EFIT  were 
tested  against  several  sets  of  hand  calculations.  These 
programs  were  found  to  be  accurate  to  within  the  number  of 
significant  digits  supported  by  single  precision  computer 
operations.  These  results  are  not  presented. 

Val idation  of  the  Program  DIATOM 

The  program  DIATOM  was  tested  against  the  analytic 
solutions  of  the  single  harmonic  oscillator.  The  harmonic 
oscillator  was  chosen  since  the  lower  vibrational  states  of 
most  diatomic  molecules  are  nearly  harmonic.  The  analytic 
energy  levels  and  wave  functions  of  the  harmonic  oscillator 
are  presented  in  Appendix  J. 

DIATOM  was  tested  using  atomic  units  where  fi  =  (i  =  1. 
The  potential  energy  model  chosen  was: 

V)  =  plr2  +  *2  1801 

The  parameter  p-^  is  related  to  the  scaling  constant  a  of 
Appendix  J  by: 

Pi  =  (ha)2  (81) 

2M- 
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and  to  (J  and  (j,  by: 


=  ipci2  (82) 

The  parameter  P2  merely  shifts  the  minimum  of  the  potential 
curve.  This  parameter  was  used  to  make  the  problem  of 
finding  the  correct  parameter  set  (p^,p2)  more  difficult  for 
MINUM . 

The  solution  of  the  wave  equation  by  DIATOM  shall  be 
referred  to  as  WAVE.  This  means  that  the  input  records  RM=0 
and  RM  =  1  were  used  in  unit  11.  The  use  of  DIATOM  to  search 
through  the  parameter  space  to  find  the  correct  parameter  set 
is  referred  to  as  MINUM.  This  means  that  unit  11  contains 
the  records  RM=1  and  RM=0.  WAVE  was  tested  first. 

The  WAVE  version  of  DIATOM  was  executed  with  p-^  =  2  and 
p2  =  0.  This  means  that  a  is  one  and  W  is  two.  Then  a 
variety  of  grid  resolutions  were  used.  Each  grid  started  at 
r  =  -6  and  ended  at  r  =  6.  The  number  of  elements  in  the 
grid  were  varied  from  five  to  ninety.  The  results  of  these 
tests  are  presented  in  Tables  IV-1,  IV-2,  and  IV-3. 

The  energy  eigenvalues  for  the  analytic  and  numeric 
(calculated  by  DIATOM)  cases  of  25,  50,  75,  and  90  element 
grids  are  presented  in  Table  IV-1.  Notice  that  the  numerical 
eigenvalues  approach  the  analytic  eigenvalues  as  the  grid 
becomes  finer.  The  solutions  of  DIATOM  converge  to  the 
analytic  solution  as  the  grid  element  size  becomes  smaller. 
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Table  IV-1 


i 


Y 


Y 


% 


Analytic  vc.  Numeric  Eigenvalues 
of  the  Single  Harmonic  Oscillator 

NUMBER  OF  GRID  ELEMENTS 


V 

ANALYTIC 

25 

50 

75 

90 

0 

1.0 

1 .000013 

1.000000 

1.000000 

1.000000 

i 

3.0 

3 .000106 

3.000002 

3.000000 

3.000000 

;* 

9.o 

9 .  ooo9  31 

5.000011 

5.00000 1 

5.000000 

i 

7.0 

7.001212 

7.000032 

7.000003 

7.000001 

4 

9.0 

9.002714 

9.000077 

9.000008 

9.000003 

5 

11  . 0 

11.00311 

11.00016 

11 . 000O2 

11.00001 

6 

13.0 

13.00931 

13.00028 

13.00003 

13.00001 

7 

15.0 

15.01265 

15.00047 

15.00005 

15.00002 

Q 

17.0 

17  *  02656 

17.00073 

17.00008 

17.00003 

9 

19.0 

19.01727 

19. 00108 

19.00012 

19.00004 

10 

21  .0 

21.07060 

21.00153 

21.00018 

21.00007 

11 

23.0 

23.01505 

23.00211 

23.00025 

23.00009 

12 

23.0 

25.12566 

25.00282 

25.00034 

25.00013 

J  3 

27.0 

27-07469 

27.00368 

27.00045 

27.00017 

14 

29 . 0 

29-11559 

29.00471 

29.00059 

29.00022 

i  'J 

31.0 

31.19916 

31.00593 

31.00075 

31  .00028 

i  6 

33.0 

33-16257 

33.00734 

33.00094 

33.00035 

.1? 

39.0 

35.22109 

35.00898 

35.00116 

35,00044 

1ft 

37.0 

37.32730 

37.01085 

37.00112 

37.00054 

19 

39 . 0 

39.37048 

39.01297 

39.00171 

39.00065 

20 

41 . 0 

41.39660 

41.01537 

41.00205 

41.00078 

21 

43.0 

43. 50121 

43.01805 

4  3 . 00  24 3 

43.00093 

‘  >  J 

4  .  .  . 

4  9 . 0 

45 • 62068 

45.02103 

45.00285 

45,00110 

23 

47.0 

47.70044 

47.02432 

47.00332 

47.00127 

24 

49.0 

49.77736 

49.02793 

49.00380 

49.00144 

Table  IV-2 


NUMBER  OF 
GRID  STEPS 


An  Overview  of  Eigenvalue  Accuracy 

NUMBER  OF  CALCULATED  EIGENVALUES 
WITHIN  X#  OF  ANALYTIC  VALUE 

CPU 

SEC.  0,005#  O.Ol#  0.05#  0.1# 


Table  IV-3 


NUMBER  OF 
GRID  STEPS 

An  Overview  of  Wave  Function  Accuracy 

%  AGREEMENT  BETWEEN  WAVE  FUNCTIONS 
v=0  v=8  v=l6 

AT  r-0 
v  =  24 

10 

0.  ?0 

* 

* 

* 

20 

0.14 

* 

* 

* 

30 

0.03 

4.40 

* 

* 

40 

r — 1 

o 

o 

1.18 

5 . 44 

* 

50 

0.003 

0.4? 

1.91 

4.82 

60 

0.002 

0.22 

0.86 

2.05 

70 

0.001 

0.12 

0.45 

1.03 

80 

0.001 

0.0? 

0.27 

0.59 

90 

0.005 

0.05 

0.17 

O.37 

*  The  -corresponding  eigenvalue  did  not 
meet  the  0.1%  accuracy  criteria. 
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This  is  expected  since  the  problem  was  constructed  for 
convergence  (17:114-115). 


Table  IV-2  is  an  overview  of  the  accuracy  obtained  for 
grids  of  different  resolutions.  DIATOM  computed  the  25 
lowest  energy  eigenvalues  to  better  than  or  equal  to  0.005% 
for  a  grid  of  90  elements  (r  =  -6  to  6).  On  the  other  hand, 
all  25  eigenvalues  were  computed  to  within  0.1%  accuracy  for 
a  lower  resolution  grid  of  45  elements  (r  =  -6  to  6).  The 
higher  resolution  grid  is  over  26  times  more  expensive 
(13  77.8  vs.  51.3  CPU  seconds)  to  use  than  the  lower 
resolution  grid.  This  is  the  trade  off  the  user  must  come  to 
terms  with.  The  computation  cost  must  be  considered  when 
specifying  the  desired  accuracy. 

Table  IV-3  presents  an  indication  of  the  accuracy  of  the 
wave  functions  calculated  by  WAVE  (DIATOM).  The  approximate 
wave  functions  vary  the  most  from  the  analytic  wave  functions 
at  the  center  of  peaks  or  troughs.  Therefore,  the  wave 
functions  are  compared  at  r  =  0,  the  middle  peak  or  trough  of 
the  even  (v  =  0,2,4,...)  wave  functions.  This  comparison  is 
only  an  indication,  not  a  true  measure  of  the  wave  function's 
accuracy.  No  value  was  given  in  Table  IV-3  if  the  accuracy 
of  the  eigenvalue  was  greater  than  0.1%,  or  if  the  wave 
function  did  not  have  the  correct  form.  Again,  notice  that 
the  approximate  solutions  calculated  by  DIATOM  converge  to 
the  analytic  solutions  as  a  finer,  more  expensive  grid  is 


used.  Also  notice  that  the  more  accurate  wave  functions  are 


the  lower  state  wave  functions.  In  all  cases  the  v  -  24  wave 


function  was  the  least  accurate.  The  best  accuracy  obtained 
for  v  -  24  at  r  =  0  was  0.37%  while  0.005%  accuracy  was 
obtained  for  v  =  0.  Better  accuracy  should  be  possible  with 
a  grid  of  more  than  90  elements. 

The  v  =  0,  3,  7  and  9  wave  functions  are  plotted  in 
Figure  IV-1.  The  dotted  line  is  the  analytic  wave  function 
and  the  solid  line  is  the  numeric  wave  function  for  a  65 
element  grid.  Each  numeric  wave  function  varied  no  more  than 
0.3%  (in  the  sense  previously  discussed)  in  absolute  value 
from  the  analytic  wave  function.  The  v  =  0  and  7  wave 
functions  differ  from  the  analytic  by  a  minus  sign.  This  is 
an  artifact  of  the  QR  method  used  by  the  IMSL  eigenvalue 
routine  EQRT2S.  Each  eigenvector  is  unique  to  within  a  minus 
sign  (11).  The  inversion  of  a  wave  function  is  of  no 
consequence  though,  since  the  Franck-Condon  factor  depends  on 
the  square  of  the  inner  product. 

Next,  DIATOM  was  tested  using  MINUM  (RM  =  1,  RVJ  =  0)  to  see 
if  the  correct  potential  energy  parameter  set  could  be  found. 
DIATOM  was  run  four  times  in  this  manner  using  the  parameters 
and  search  limits  of  Table  IV-4.  After  the  fourth  run,  MINUM 
chose  the  parameters  p-^  =  2.0004  and  p2  =  -0.002.  These 
values  are  close  to  the  correct  values  of  p^  =  2.0  and  p2  = 
0.  A  closer  fit  is  possible  by  running  MINUM  (DIATOM)  again 
with  a  narrower  search  region  and  more  grid  elements. 
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-6  -4  -2  0  2  4  6 

V  »  0 

Fig.  IV-la.  Selected  Analytic  and  Numeric  Wave 

Functions  of  the  Single  Harmonic  Oscillator 
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-6  -4  -2  0  2  4  6 

v  -  7 

Fig.  IV-lc.  Selected  Analytic  and  Numeric  Wave 

Functions  of  the  Single  Harmonic  Oscillator 
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Table  IV-4 


RUN 

NUMBER 

1 


2 


3 

4 


Test  Runs  of  DIATOM  Using  1WINUM 

P2-0.0 


Correct  Set  P-^-2.0 


KEY  WORDED 
RECORDS 

NE;  20 
IS "100 
Pl:-2 .  6 
LI -0.1 
Ui  o.o 
P 2-0.0 
L2--10.0 
U2=10.0 


RESULTS 


Chosen  seti  P^-l*98o3 


P2=0.0221 


.urn 


of  residuals  squared  =  1.6  x  10 


*  Pl=1.98 
LI =1.90 
U1--2.10 
P  2-0 . 02 
L1--0 . 02 
Ul=0.03 

*  NE-30 


Chosen  seti 


*  NE=40 
IS  =  60 
P I  2 . OOn 

LI =1.99 
UI =2 . 01 
1  2=-0.01 
L2--0.02 
U2=0 . 02 


P. -1.9800 
P  2=0 . 0200 


Sum  of  residuals  squared  =  ?.8  x  10 

Chosen  sett  P^-2.0061 
P2=-0.0196 

Sum  of  residuals  squared  =  6.1  x  10 
Chosen  sett  P^=2.0004 

P2  -o.uoiv 

Sum  of  residuals  squared  =  8.8  x  10 


-4 


-6 


*  All  other  parameters  the  same  as  the  previous 
run. 
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DIATOM  appears  to  be  capable  of  solving  the  wave 
equation  accurately.  The  level  of  accuracy  depends  on  the 
grid  used.  Also,  DI ATOM  seems  to  be  able  to  find  the  best 
set  of  potential  energy  parameters  for  the  system  of 
interest . 

DIATOM  is  best  used  in  a  two  step  process.  First,  low 
resolution  grids  and  large  numbers  of  steps  through  M1NUM  are 
used  to  narrow  the  search  region  for  each  parameter.  Then 
the  resolution  of  the  grid  is  increased  while  parameter 
limits  are  reduced  until  an  acceptable  fit  of  energy  levels 
is  achieved.  Then  MINUM  is  turned  off  and  the  number  of  grid 
elements  increased  so  DIATOM  can  calculate  accurate  wave 
f  unct ions . 

Va 1 idation  of  the  Program  FCFACT 

The  program  FCFACT  was  tested  only  on  the  wave  functions 
of  the  same  harmonic  oscillator  problem.  True  Franck-Condon 
factors  were  not  calculated  since  the  wave  functions  used  all 
came  from  the  same  potential.  Therefore  the  factors  computed 
by  FCFACT  were  the  square  of  the  inner  product  (i  =  0 
to  24;  j  =  0  to  24).  For  true  Franck-Condon  factors,  arid 
i|tj  would  belong  to  different  potentials.  Also,  each  set  of 
wave  functions  would  be  calculated  by  runs  of  DIATOM. 

Tire  square  of  the  inner  product  was  as  expected  since 
the  harmonic  oscillator  wave  functions  are  orthogonal.  If  i 
=  j  then  the  inner  product  is  1.0.  Otherwise  it  is  zero. 


These  results  were  achieved  with  a  grid  of  sixty  elements. 
The  square  of  the  inner  product  was  non-zero  (.001  to  .007) 


between  high  vibrational  wave  functions  (v  >  15)  for  a  45 
element  grid.  These  non-zero  values  occured  only  between 
and  ij,  j  when  i  =  j  _+  2. 

The  square  of  the  inner  product  was  also  computed  using 
Simpson's  integration  rule  over  1001  points  for  v  =  0  to  9. 
The  spline  coefficients  were  used  to  interpolate  the  1001 
points  along  the  grid.  All  Simpson's  rule  values  agreed  with 
those  of  FCFACT. 

The  non-zero  values  arose  since  the  numerical  wave 
function  only  approximates  the  analytic  wave  function.  The 
approximation  could  be  made  worse  if  a  bad  cubic  spline  fit 
is  made  to  the  numerical  wave  function.  Still  though,  FCFACT 
gives  good  results  for  high  resolution  grids  of  60  elements 
or  more. 

Program  DIATOM  Benchmark 

Many  runs  of  DIATOM  were  made  in  order  to  characterize 
its  CPU  time  requirements.  These  runs  are  summarized  in 
Table  I V- 5 .  The  typical  CPU  requirements  of  the  harmonic 
oscillator  are  presented.  The  first  column  indicates  how 
many  grid  elements  were  used.  Each  grid  began  at  r  =  -6  and 
ended  at  r  =  6.  The  second  column  indicates  how  long  the 
program  took  to  solve  the  wave  equation  (RM=0,  RW=1).  MINUM 
was  not  allowed  to  run  for  this  data. 
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■irfwiidt 


The  third  column 


RD-A151  765 
UNCLASSIFIED 


COMPUTER  MODELING  OF  VIBRATIONAL  ENERGV  LEVELS  OF  2/2 

POTENTIAL  LASER  CANDIDA.  .  <U>  AIR  FORCE  INST  OF  TECH 
MRIGHT-PATTERSON  RFB  OH  SCHOOL  OF  ENGI.  .  P  H  OSTDIEK 
DEC  84  AFIT/GEP/PH/84D-6  F/G  28/8  NL 


Table  IV-5 


An  Overview  of  Program  DIATOM 'a  CPU  Uce 


- Rlv|  ”l ,  HW-O,  13-1 - 

NUMBER  OF  --[(M  O,  RW-  l  --  TOTAL  ONE  PASS  OVERHEAD 


GRID  ELEMENTS 

TOTAL  CPU  SEC 

CPU  SEC 

CPU  SEC 

CPU  SEC 

5 

1.9 

1.83 

0.11 

1.72 

10 

3.3 

2.64 

0.31 

2.33 

15 

5.2 

3.99 

o.e>4 

3-35 

20 

9.0 

8.18 

1.14 

7.04 

25 

14 . 2 

12.26 

1.82 

10.44 

30 

17.0 

12.40 

2.75 

9.65 

35 

25.0 

17.05 

3.98 

13.07 

40 

35-0 

24 . 54 

5.84 

18.?0 

45 

51.3 

45,86 

11.15 

34.72 

50 

76.2 

IbO .  64 

39.83 

120.81 

55 

12 1  .1 

369.35 

92.05 

277.30 

60 

180,4 

994.48 

180. 03 

814.45 
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contains  the  CPU  time  required  for  DIATOM  to  run  when  MINUM 
was  allowed  only  one  step  (RM=1,  RW=0,  IS=1).  The  fourth 
column  contains  the  CPU  time  between  calls  of  the  function 
FUN  by  the  subroutine  MINUM.  This  is  the  approximate  time 
for  MINUM  to  take  one  step,  and  is  referred  to  as  the  "one- 
pass"  time.  The  fifth  column  contains  the  remainder  of  the 
CPU  time  for  DIATOM  to  run.  This  is  the  overhead  associated 
with  the  rest  of  the  program  and  is  referred  to  as  the 
"overhead"  time. 

To  approximate  how  long  a  run  of  DIATOM  will  take  (RM=1, 
RW=0)  use: 

CPU  time  =  (overhead  time)  +  IS (one-pass  time)  (83) 

where  IS  is  the  number  of  steps  MINUM  is  allowed  to  take.  Eq 
(83)  calculated  the  CPU  time  to  within  30  seconds  for  the 
runs  of  DIATOM  discussed  previously. 


V. 


Summary  and  Recommendations 


Summary 

Four  programs  were  written  to  be  used  as  a  set  to 
calculate  the  Franck-Condon  factors  between  two  electronic 
states  of  a  diatomic  molecule.  The  factors  calculated  are 
for  v  =  0  to  24.  These  programs  are  DUNHAM,  EFIT,  DIATOM, 
and  FCFACT.  DUNHAM  calculates  approximate  energy  levels 
using  the  Dunham  equation  {Eq  (1))  and  coefficients.  Program 
EFIT  uses  a  least  squares  technique  to  find  a  set  of  energy 
levels  represented  by  spectroscopic  data.  Both  programs  are 
used  to  build  the  input  file  for  DIATOM.  This  file  contains 
the  observed  energy  levels  to  which  the  MINUM  portion  of 
DIATOM  will  fit  a  potential  energy  curve.  DIATOM  uses  a 
finite  element  technique  to  solve  the  wave  equation  and 
calculate  wave  function  values.  These  wave  function  values 
are  used  by  FCFACT  to  calculate  Franck-Condon  factors. 
DIATOM  has  derived  very  accurate  wave  functions  for  the 
single  harmonic  oscillator.  Also,  FCFACT  has  calculated  the 
proper  Franck-Condon  factors  between  these  wave  functions. 
This  program  set  promises  to  be  an  inexpensive  and  accurate 
alternative  to  the  RKR-IPA  program  set.  However,  it  still 
needs  testing  and  modification. 


Recommendations 


1.  Rewrite  DIATOM  to  solve  the  wave  equation  in 
dimensionless  form.  This  reduces  the  chance  of  mathematical 
operations  exceeding  storage  limits  (underf low/ overf low 
errors)  of  the  computer.  Also  by  adding  a  small  unit 
conversion  subroutine,  the  user  could  use  whatever  units 
desired. 

2.  The  execution  time  of  DIATOM  can  be  reduced  by 
finding  faster  eigenvalue  routines.  Specifically,  Dr. 
Shankland  has  written  one  routine  which  should  be 
investigated,  and  used  if  faster  than  IMSL's  EQRT1S. 

3.  DIATOM  and  FCFACT  need  to  be  tested  against  the 
analytic  solutions  of  the  Morse  potential  wave  equation. 
This  should  be  done  as  this  work  did  for  the  single  harmonic 
oscillator. 

4.  Compare  the  results  of  this  program  set  and  the  RKR- 
I P A  set  for  lead-oxide  and  lithium  hydride.  The  RKR-IPA 
results  are  available  in  Pow's  work  (16). 

5.  Move  the  program  set  to  the  DEC  VAX  11/780  and  run 
under  the  UNIX  operating  system.  This  version  of  the 
programs  would  be  more  transportable  since  both  DEC  equipment 
and  the  UNIX  operating  system  are  popular  in  laboratories. 
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Program  DUNHAM  Flow 


Program  DUNHAM 


Program:  DUNHAM 

Version:  84  .  05  .  i  3 

Author:  Paul  H.  Oililu-l; 

Air  Pure  -  institute-  pi  T^tiii^ioj/ 

Ur  ttr  tur.  An  Force  L'uit- ,  Oh 

Dtficr  .ptior.:  inis  ,ji  Oji'cin  calculates  approximate  energy  levels 

tor  diatomic  molecules  using  the  Dunham  equation. 

I/O  logical  unit  6  --  output  listing  tile 
11  -  -  i nput  file 

13  --  output  listing  file 
( Same  as  ui.it  o) 

1  <1  --  energy  level  output  file 


PROGRAM  main 


CHAftACTER*2  NUMDER <0: 25) 

CHARACTER*30  NEnDER 

INTEGER  1,  J,  LVLLIMT,  VID,  VIDLMT,  ROT,  KOTLMT ,  EVlii(67o), 

EPOT  (  67’o ) 


REAL 


SUM,  T < G: 25, 0: 25) ,  Y (0:9, 0:9),  RK01  (O: 25)  , 
RV In ( 0 : 25 )  ,  ELVL  v  o76 )  ,  DEQUIL 


COMMON  /HD  r.  /  HEADER 

COMMON  /DAT/  Y,  DEQUIL,  VIELMT,  ROTLi'lT,  LVLLMT 


DATm  NUrlOER  /  ’  00  ’  , 

’  Oi  ■ 

,  ’02’ 

,  ’  03’ 

,  ’  04 ’  ,  ’ 05’  ,  ’ Go’  , 

’  10’  , 

’ll' 

,’12” 

,  ’  13’ 

,  *  14*  , ’ 15’  ,  ’ 16’  , 

’  20  ’  , 

'  21  ’ 

t 

*  O’?  » 

, 

,  ’24’  ,  >25’ / 

07 ’ , ’ 08 ’ , ' 09 ’ , 
17’ , • to • , ’ iy ’ , 


DATA  RRGT 


/  0.0,  1.0,  2.0,  3.0,  4.0,  5.0,  6.0,  7.0,  8.0,  9.0, 

10. 0,11. 0,12.u,13. 0,14. 0,1 5. 0,16. 0,17. 0,18. 0,19.0, 
20.0,21.0,22.0,23.0,24.0,25.0/ 


r>M r a  r-:vifa 


/  0.0,  1.0,  2.0,  3.0,  4.0,  5.0,  6.0,  7.0,  8.0,  9.0, 

10.0,11.0,12.0,13.0,14.0,15.0,16.0,17.0,18.0,19.0, 
20.0,21.0,22.0,23.0,24.0,25.0/ 


Open  the  output  listing  file,  print 
the  header  ,  and  star  t  t  uri  stats. 


CALL  HDRDLIM 


-Re&d  the  input  tile  (unit  II  J  ter 
the  contr  ol  pm  uinctcf  a  alia  date. 


C Mi.  L  RDRDUN 


WRITE  (13,1 303 ) 

FORMAT  (  '  1  tiler  l  es  T  (  v  ,  J  > 


(T-0  it  potential  in  i  n  linum)  ’  ,  /  ) 


-Find  tlie  Dunnmn  Equation  energy 
TiVIB.ROT)  where: 

VIB  is  the  vibrational  quantum 
number  of  the  current  level 
ROT  is  the  rotational  quantum 

number  or  the  current  level 


DU  *40  0 IE-0 ,01  BLI'IT 
I/O  30  ROT  =  0 ,  R0TLMT 
SUM  =  0 


■The  sum  (do  loop)  indicies  may  vary 
from  0  to  9.  This  yields  the  first 
100  terms  of  the  Dunham  Equation. 


DO  30  1-0,9 
DO  10  1-0.9 
IF  i 1  . EG 


IF  ii  .EG.  0  . mHD.  J 
SUM  -  SUM  t  V  l  I  ,  J  ) 
El  "C 

IF  (I  .EG.  0  .AMD. 


.EG.  0)  THEN 


IF  (I  .EG.  0  .AMD.  J  .NE.  0)  THEN 

SUM  -  SUMt ( Y ( I , J ) KKkOT (ROT) **J* (RR0T (ROT ) +1 . 0) **J ) 


IF  (1  . HE .  O 


J  .  EU.  0)  THEN 


=  SUM  +  (Y(I,J>  *  (RVIBCVIB)  +  G.S)**I) 

—  wUi’H  (  V  (  1  ,  J  )  it  (  HV  IB(VIB)  »0.  5)  ttlfl 
arKOT  (ROT)  **J* (RROT  (ROT)  +1 . 0)  **J  ) 


END  IF 
END  IF 
END  IF 
CONTINUE 
CONTINUE 


-Transfer  the  accumulated  sum  to  T 


T  (GIB,  ROT  i  -  U UN 

WRITE  (13,130*4)  VIE,  ROT,  T  (VXD, ROD 
FORMAT  (’  T  (’  ,13,  ’  ,  ’  ,13,  ’  )  —  ’  ,  G1 1> .  7 ) 
CONTINUE 
CONTINUE 


Transfer  the  calculated  energy 
levels  T  i U I B , R 0 T /  to  ELVL ( I )  . 


I  -  0 

DO  aO  VIL-Oj  VltsLMT 
Du  DO  ROT-O,  ROTLMT 
i  -  I  ■*  1 

ULVL  (  X  »  -  1  (VXD,  KG  i'  ) 

EV I 6 ( I >  =  V 1  i. 

EKG1 ( I )  -  ROT 
DO  COM  riNGt 

oO  CONTINUE 

C  - — - — . . — - -  - - oh  1  t  t  tiler  e  h  e  r  0  /  ltrv/dlo  b»Q  that  t  n  e 

C  lowest  tVID~0  de  ROT^U)  lb  at 

C  O  -  (the  dissociation  energy) 

Du  ~’0  I  -  i  ,  L  OLi-h.  I 

LLVL  (  I  )  --  :■  LVL  (  1  )  -  DiIGtU  I  L 

70  C  uWT  I  i'itlc 


Ui-  i  IE  C",i:01)  HEADER 

1  30  1  FGM’IA'i  l  *  1  Dunham  Equation  tnei  y/  Levtlb  ’  ,hJG,/) 

DO  80  I-l.LVLLMT 

Ui->  I  TL  (ID,  IV  OX:  >  LVIDCI),  EkU)  (  I  )  ,  LLVLil) 
l-’OL  rU.-Tini  (  ’  w-’.i-i,’  J-*  ,  14,  IX, G15.?) 

:o  LuNiiNUe 

C  - - -  - - - Hut  the  energy  levels*  into  an 

C  Output  tile  (unit  14). 


'•REN  iullll-1-O 
DO  9u  I  -  1  ,  LV LLI'I  I 

Ur-  'I  TE  (i -1,1401)  Null'll:  i£l<  (1*1),  ELVL  (  I  ) 
i  -i  -0  1  FORMA  T  (  ’  ;-vr  ,  A2,  ,Glt».7,  ’  5  ') 

COM TIN JE 

CLOCu  (uHIT-14) 


C  -  -  -  - - C  1  osie  the  output  Tile  and  run 

C  btjtlbtlCS. 


CALL  TKLDUN 

9999  END 
♦ADD, HDRDUN 
♦ADD , RDRDUN 
■4-mDD,  TI-uDUN 
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C  frjji'oa:  SUBROUTINE  tiDKDUN 

C 

C  Vt rrizian:  C*l.  OB.  13 

C  Author :  Paul  H.  ust-dieL 

C 

C  Air  Force  Institute  ot  Technology 

C  Wr i ght-Patter son  Air  Force  Base,  uH 

0 

C  Description:  HJJRDUN  opens  the  output  listing  tile  (logical  units 

C  13  arid  a)  and  Starts  CPU  and  wall  time  use 

C  statistics. 

C 

C - - - 

I-UfROUTIHE  HDRDUN 

CHAR  AC  T  £P  VER'i'N 

CHAR  AC  TL‘  Kit  1  3  PCM 

1  HTLot.K  *3  i  DnTE  ( 3  >  >  1IXNE(3),  1USLK(4) 

t’LK-l'l  =  ’  o-l.0c.13’ 

f  CM  -  ’  Citr/84D-o/  1.2* 

C  - Initiate  1 1 1 e  run  statistics. 

CALL  ETINE 
CALL  ST  I ME 


C  - Get  the  current  date,  time,  and 

C  user  name  tor  output  on  the  header 


CnLL  DATEilliAiffc) 

call  riMEdTiMfc) 

C  ALL.  oCLKNG  (  I UsER) 

C  - - Open  the  output  listing  tile  and 

C  write  out  the  header 

C  rLi  l  (  Un  IT-  i  3  ) 

OPEN  (UNI  1-6) 

UPirt  (13,1301)  IUGER,  VERCN,  IDATC,  PCN,  I T 1 ML 
1..1  FORMAT  (’1  User:  ’ , A A3 , Tb l , ’ A i r  Force  Institute  ot  Technology’, 

i  T  1  1  0  ,  '  (/ e r  s  i  o n  :  '  ,  AC, 

*  Date:  ’  ,3A3,Tll-4,  ’PCN:  ’,A13, 

>  Til..-:  ’  ,  3A3  ,  lot ,  ’  DUNHAM  ’,/,/) 

VWI  )AL  f  Ur.  1 1 

L  .'l  l) 
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ENDIF 

CONTINUE 

ELSE 

IF  <SC0RE(2:5)  .EG.  ’ LULS ' )  THEN 

- Find  the  correct  state  number  xx  at 

LVLS.^x  where  xx  is  ’01’  to  '10'  and 
corresponaa  ta  STATxx  above. 

DO  250  1=1,10 

IF  (SC0REC6:/)  .Eu.  DIOITS(I))  THEN 

- Count  how  many  levels  are  associated 

with  the  state  xx  ana  record  them  in 
STATE ( xx , count ) 

z>  0  0  L  =  V 
UU  2AG  j  =  i  ,  lo 

DO  2-0  L=CC0L,7i 

IF  (SCORE (L:L)  .WE.  ’  ’)  THEN 

IF  (WINS! (I)  .LT.  25)  THEN 
ITiNST(I)  =  NihST(I)  t  1 
READ  (SCORE (L: (L+l)),’ (12) ’/ 

ST  A  f  E ( 1 , N I NS  f ( 1 )  ) 

SCOL  =  l-  v  2 
ENDIF 
ENDIF 
Co.il  I NUE 
CON  I  I  HUE 
ou  ro  so 
cl  (D  I  F 
CuN  i  INUlI 

. F  (SCORE (2: 7)  . EG .  ’SHIFTS’)  THEN 

- SHIFT  » a  the  vaxue  deaired  tor  the 

loweat  level  ot  the  loweat  atate. 

Kt.AU  (  -■  C  UKt ( V : 2a )  ,  ’  (  E  1 5 . 7  )  ’  )  SHIF  1 
EL  ^2 

- - Look  tor  a  keyword  ABBCDD  where: 

A  E  C  match  valuee  in  STLBL 
BB  E  DD  match  valuee  in  STATE. 

FOUND  =  .FALSE. 

- - Look  tor  part  A. 

Du  S1U  1-1,10 

IF  '-.CORE  (2:  2)  .EG.  STLbL ( I ) >  THEN 

CUK-  i  ( l :  i )  =  scorl.  (2:2) 

I  Li  t  =  I 

Uvli.i 

G  ei  i  i  1 1  )OL 
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DO  0  Id , 10 
STLBL  l  I  )  = 

NINSTd)  =- 
DO  2  J-  1  , 2 1- 

OTA  TO  (  I  ,  J  )  =-  O 

COM  1’  I  NUE 
C jNTIWUE 

L'Ci  5  I  -  1 , 200 
DO  -4  J= 1,250 
L  I  ME ( I , J  )  =  O 

WEIGHT  (1,11)  -  0 


-t  CONTINUE 

5  CONTINUE 

C  - up u  n  the  input  data  Tile. 


OPEN  ( UN I  T=  1 1 ) 

WRITE  (13,1001) 

1301  FORMAT  </,’ O’ , T31 ,’ Input  Data’,/,’  ’,T31,’ -  - ’) 

DU  00  K- 1,10000 

C  - Transfer  a  record  trum  the  input 

C  Tile  to  the  buTfer  SCORE. 

READ  <  1  1  ,  1  101  ,  END-81  )  SCORE 
1101  FORMAT  ( A 72 ) 

WRITE  (13,1302)  SCORE 

1302  FORMAT  ('  ’  ,  A’’2) 


c  - ]r  SCORE  (111)  te  a  ’  then  SCORE 

C  should  contain  data  and  a  valid 

C  le/word.  So,  see  i T  SCORE (217) 

C  lues  contain  a  valid  keyword.  IT  it 

C  does.,  t r  an sTer  the  data  Trom  SCORE 

C  to  the  input  variable. 


IF  (SCORE* 111)  .EQ.  ’>’>  THEN 

IF  (SCORE  (2:5)  .EC.  ' STAT’)  THEN 

3  - Fit  id  the  correct  state  number  STATxx 

C  where  xx  ia  ’01’  to  ’10’. 

DO  22u  1-1,10 

IF  ( SCORE (6:7)  .EQ.  DIGITS(I))  THEN 
DO  210  J  =  9 , 72 

IF  (SCORE ( J : J  >  . NE.  *  ’)  THEN 

STEBL ( I ) (1:1)  =  SCOREIJIJ) 

SETCNT  =  SETCNT  ♦  1 

GO  TO  60 
END  I F 

210  CONTINUE 
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SUBROUTINE  RDREFT 


VcfUloiw  Sd  .  OU  .  13 

C  Author:  Paul  H.  OLtdiul; 

C 

C  Air  Force  Institute  of  Tecluiology 

C  Wr  i  jlit-Pattur  &uri  Air  Faroe  Base,  QH 

C 

C  or  1  p  t  j  or. :  Tliit  routine  opens  an  input  tile  (unit  11)  and  reads 

all  r«--carda  within.  Each  record  read  is  written  to 
•l  the  output  listing  (unit  13).  Bata  records  are 

inerted  by  «  ’>’  in  column  1.  These  records  contain 
d.-ta  referenced  by  a  single  hey  word.  All  other 
C  records  are  considered  comments.  This  routine  uses 

S  ah  internal  read,  eg  READ  (SCORE (5: 19) , ’ (Elb. 7) ’ )  X 

reads  from  columns  5  to  19  of  the  character 
variable  SCORE  using  the  edit  descriptor  Elb. 7  into 
C  the  real  variable  X. 


SUBROUTINE  RDREFT 

CHARACTER* 1  CURST,  STEEL (10) 

CHARACTER*2  NUMB 2 6 ( 26 ) ,  DIGITS (10) 

CHARACTER*/’  P  UkM 
CHrtRAC  IEK*.'S  .,CO:£ 

INTEGER  I,  LI,  J,  L2,  K,  SIZE,  FCOL,  LCOL,  NINSKIO), 
CURLVL,  LVLNG,  ILBL,  SETCNT ,  STAT£(10,2b) 

cOGICAL  FOUND 

REAL  L INE <2b0, 2SC) ,  WE I GHT l 250 , 250 ) ,  SHIFT 

summon  /indata/  line,  weight,  size,  shift 

COMMON  / LABELS /  a ILBL ,  CURST 

COMMON  / LEVEL  0 /  NXNST,  STATE,  CURLVL,  SETCNT 


2Alr>  DIGITS 

/  ’  01  ’ 

,  ’02' 

,  ’03' 

, ’Od’ , 

’Ob’, 

'On', 

■0.”  , 

’OS’ , ’09’ , 

’  10’  / 

DATA  i IUMB26 

/  ’  00’ 

,  ’01' 

1  ,  ’ 02 ’ 

, ’03’ , 

’  o-t  ’  , 

’Ob’, 

’Ob’, 

’0,”  ,  ’ OS’  , 

’09’  , 

t- 

’  10  ' 

,’11’ 

,  '  12  ’ 

, ’ 13’ , 

’  1R’  , 

lie. 

lu  , 

’  16’  , 

’  17'  ,  ’  Id  , 

'  19’  , 

♦ 

1  20  ’ 

f  *-)  |  1 

J  * 

f 

»  | 

»  9 

’2A’  , 

’  2b  ’  / 

E 


Initialize  the  input  Variables. 


s  SIZE  --  0 

SET CN  r  -  0 
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c - 

c 

C  Pruji- LUbkUU  7  INE  HJJKliFT 


C  Version:  .  Oii  .  13 

C 

C  Author  :  Paul  H.  OstdieL 

c 

C  Air  Force  Institute  of  I  ec  huo  1  oy  y 

C  Wr  l  yh  t.  -  Pat  t  er son  Air  Force  Base,  uH 

C 

C  Description:  FiUKLI-' 1  opens  true  output  i  lbtiinj  tile  (logical  unite* 

C  1j  ana  oJ  end  starts  CPU  and  wall  time  use 

C  a  tat  i  s  1 1  c  s . 


CUERuUTINE  HDREFT 

CHAKAC  I  ER*U  VERSN 
CHARACf EK*13  PCN 

INTEGEF<*3  IDATE  (  3  )  ,  1TIHE(3>,  IUSER(A) 


VERSN  •-=  ’U-l.0w.13’ 

PCN  -•  ’  GEP/UAD-6/  1  .  1  ’ 

C  - Initiate  the  run  statistics*. 


CntL  LT  X  NE 
Call  STINE 

C  - - - Get  tire  current  date,  tune,  and 

C  user  name  for  output  on  the  header 

CALL  DA it (I  DATE) 

CALL  1  INC.  <  ITIME) 

CALL  USEKNO ( I USER ) 

C  - Open  t n e  output  listing  file  ana 

C  write  out  the  header 

open  (uun-13) 

OPEN  ( UNI T=6 ) 

WRITE  (13,1301)  XUOEK,  VERSN,  I  DATE,  PCN,  ITIME 
1  ;  1  i  I-  UK  MAT  (’I  User:  ’  ,  -1A3,  T51 ,  ’  Ai  r  Furce  Institute  of  Technology’, 

+  T110, ’Version:  ’,A0, 

«  /,’  Date:  ’ , 3a3 , TI 14, ’ PCN:  ’,A13, 

*  /,’  Tunc-:  ’  ,3A3,  162,  ’ENERGY  FIT’,/,/) 

9999  RETURN 
END 
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DO  190  1=1, SIZE 

X ( I )  -  X ( I )  +  SHIP  T 

ivO  CONTINUE 

WRITE  (13, ISOS)  SHIFT 

ISOS  FORMAT  (  ’  lEner^y  levt-la  shitted  =,0  lowest  level  is  ’,015.7) 
WRITE  (13,1302)  ( X ( I ) , I =1 , SIZE ) 


C  - Write  the  energy  levels  to  the 

C  output  listing  (unit  13)  end  an 

C  output  tile  (unit  14). 


OPEN  ( UNI T— 1 4 ) 

OFFSET  -  0 
WRITE  (13,1312) 

1C  12  FORMAT  (’IThe  energy  levels  -for  each  state  are:’) 

DO  400  1=1, SETCNT 

WRITE  (13,1310)  STLBL ( I ) 

1,10  FORMAT  (’OLevels  01  state:  ’,Al) 

DO  410  J=1 , MINST ( I ) 

WRITE  (13,1311)  STATE (1,3),  X(GFFSET+J) 

1311  FORMAT  (’  ’,16,015.7) 

WRITE  (14,1401)  NUMBER (STATE ( I , J >) ,  X(0FFSET  +  J> 


1401  FORMAT  ( ’ >0 » , A2 , ’ “ * , 015 , 7 , »  i  ’) 

410  CONTINUE 

OFFSET  =  OFFSET  +  NINST(I) 
a 00  CONTINUE 

CLOSE  ( UNIT  =  14) 

C  - - Close  the  output  tile  and  run 

C  stat 1 st 1 c  s . 


CALL  TRLEFT 

9999  END 

trtDD , RDREFT 
♦ADD , LULEFT 
♦ADD, HDREFT 
♦ADD, TRLEFT 


IF  (I  . EG .  J)  THEN 

A  (  I ,  J  )  =  SORT ( A  ( I ,  J  )  -  SUM) 

ELSE 

A  (  I  ,  J  )  =  (A(I,J>  -  SUM)  /  A  (  J  ,  J  ) 

END  IF 
CON  riNUE 
CONTINUE 


-Solve  the  equation: 

LY  =  6 

where  Y  -  L  ( tr  arispose )  x  X 
and  is  stored  in  X. 


DC  l-4u  1  =  1,  (SI  ZE-  1  ) 

SUM  =  0 

DO  130  K=1 , ( 1  1 > 

SUM  =  SUM  +  ( A ( I , K )  *  a(K)> 

CONTINUE 

X(I)  =  (Bil)  -  SUM)  /  A ( 1 , I ) 
CONTINUE 
X (SILL)  =  O 


-This  is  a  consistency  check. 


DO  ISO  k=l , (SIZE- 1 ) 

SUM  •=  SUM  <■  (A  (SIZE,  K)  *  YAK)) 

CONTINUE 
WRITE  (13,1 304 ) 

FORMAT  ( ’  1  Cons i stency  check...') 

WRITE  (13,1301)  SUM,  SIZE,  B(SIZE) 

FORMAT  (’  SUM  =  ’,615.7,/,’  B(’  ,12,’)  =  ’,015.7) 


-Recover  X  From  Y  (stored  in  X) 


DO  170  I •=  (SIZE-1 ),  1, -1 
SUM  =  O 

DO  160  K= ( 1+1 ) , (SIZE- 1 ) , 1 

SUM  =  SUM  +  ( A ( K , I )  *  X(K)> 

CONTINUE 

:;  ( I )  =  (X  ( i  >  -  sum)  /  m  ( i ,  i ) 
CONTINUE 


-Write  out  X 


Wi  I  i  Z  i  1  3 ,  i  305 ) 

I  ORMA T  </,’  Unshifted  energy  levels  are...’) 
WRITE  (13,1302)  ( X ( I ) , I = 1 , SIZE > 

FORMAT  (’>’,015.7) 


■Shift  the  energy  values. 


DO  130  I =SIZE ,1,-1 
X  (  I  )  =  >:  ( I  )  -  X  (  1  ) 

CONTINUE 


T 


C  - Build  A’s  d  i  a^udal  s  o+  AX=B. 

DO  30  K=1,SI2E 
CUMIN  =  0 
DO  10  I-l, SIZE 

CUMIN  -  CUMIN  ♦  WE  1 UHT  (  1 ,  K  ) 

10  CONTINUE 

RUNOUT  -  0 
DO  20  I -^1,  SIZE 

CUMOUT  -  CUMOUT  +  WEIGHT (K, I ) 

20  C  url  r  I NUE 

A ( K , K )  =  CUMIN  +  CUMOUT 
.0  CONTINUE 

C  - build  riDn-d  i agonal  a  or  A. 

DO  CO  K= 1 >  C 1 2E 
DO  50  J =1 , (K-l  ) 

A  ( K  ,  J  )  =  -1.0  *  (WEIGHT  (K,  J  )  +  WEIGHT  (  J  ,  K)  ) 

A (J , K)  -  A  ( K ,  J  > 


50  CONTINUE 

^0  CONTINUE 

C  - Bu  lid  B 


DO  90  K-l , CI2E 


CUMOUT  =  O 
DO  70  J - 1 , S I 2E 

CUMOUT  =  CUMOUT  ♦  (WEIGHT (K,Ji  *  LINE<K,J>> 

70  CONTINUE 

CUMIN  =  0 
DO  30"-J=l  ,  SIZE 

CUMIN  =  CUMIN  +  (WEIGHT  (J,K)  *  LINE(J,K>> 

oG  CONTINUE 

B  ( K )  =  C UI'IOU  I  -  CUMIN 

90  CONTINUE 

C  - - -Solve  for  X 

C 

C  Firbt  ubr  Cholesky  decomposition  to 

C  •get  L  (stored  in  A)  such  that: 

C  A  =  L  x  L ( transpose) . 


DO  120  J-=  1  ,  S I 2E 
DO  110  I=J,CI2E 
CUM  -  0 

DO  100  K-l , ( J-l ) 

IF  (1  .  £Q .  J )  THEN 

SUM  -  CUM  +  (A(I,K>**2) 

ELSE 

SUM  =•  CUM  *  ( A ( 1 , K )  *  A  <  J  ,  K )  ) 

END  IF 

lOO  CONTINUE 
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Appendix  D 
Program  EFIT 


Program:  UK IT 

Verbiun:  04  .  Oti .  13 

Author:  Paul  H. 

Air  Force  Institute  of  I  echi.ology 
wr •  i  gh t-P«t ter son  Air  Force  base,  OH 

Description:  This  program  uses  a  least  squares  technique  to  Find 

the  set  of  energy  levels  that  best  Fit  a  set  oF 
spectroscopic  transition  lines. 

I/O  logical  unit  6  --  output  listing  File 
i 1  - -  i npu  t  File 
li  --  output  listing  File 
( same  as  unit  a) 

14  --  energy  level  output  File 


PROuRAM  MAIN 

CHARACTER*  1  CURST,  STLBLIIO) 

CHARACr£R*2  HUMBER  v  0 :  25  ) 

INTEGER  I,  J,  K,  S12E,  NINST(IO),  STATE (10, 2b) ,  CURLVL, 

*  SETCNT,  OFFSET 

REAL  A ( 2S0 , 250 ) ,  X<250>,  B<250>,  SUM,  SHIFT,  LINE ( 2 SO, 250 )  , 

♦  WEIGHT <250, 250) ,  3UMIN,  SUP10UT 

COMMON  A 

COMMON  / 1 MDAT  A /  LINE,  WEIGHT,  SIZE,  SHIFT 
COMMON  /LABELS/  STLBL,  CURST 

COMMON  /LEVELS/  NINST ,  STATE,  CURLVL,  SETCNT 


/  ’  00’  , 

’  01  ’  , 

’02’  , 

'OS’ 

,  ’  04  ’ 

, ’05’ , 

’06’  , 

’07’  , 

'OS’  , 

’  09  ’  , 

•  10’  , 

•11’ , 

’  12’  , 

’  13’ 

,  ’  14  ’ 

, ’ 15’ , ’ 

’  16’  , 

’  17’  , 

’  16’  , 

'  19’  , 

’20’  , 

’  21  ’  , 

»  r? * 

’  23’ 

,  ’  24  ’ 

, ’ 25’ / 

-Open  the  output  listing  File,  print 
the  header,  and  start  run  stats. 


CALL  HUREFT 


-Read  in  measured  line  values  and 
control  parameters. 


ChLL  RUREFT 
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Program:  SUBROUTINE  TNLDUN 

Ver a  ion:  oh . Ou. 13 

hUthor!  Paul  H.  uutdluk 

Air  Force  Iriotitnie  of  technology 
Wr  :yht  Puiterbon  Air  Furct  Lai.e,  OH 

Eo-jCi  iption:  TKLDUN  cloi,ea  t he  output  lifting  tilt?  (logical 

units  10  and  a)  and  stops  the-  CPU  arid  wall  time 
use  statistics. 


SUBROUTINE  TRLDUN 

c 

CALL  ETlt'lE 

CALL  WT I ME 

-Shut  down 

the  run  statistics 

c 

CLOSE  (UNIT- 13, STATUS- 1  KEEP’  ) 
CLOSE  ( UN IT -6  > 

-Close  the 

1 i st i ng  outputs. 

rp.ygo 

RETURN 

<  E 1 5 . 7 )  ’  )  V  ( 1 1 ,  J  J  ) 


READ  (SCORE  <6: 20) , ’ 

EMDIF 
END  IF 
EMDIF 
END  IF 
END  IF 
END  I  F 
END  IF 

30  CONI INUE 

C  - Close  the  input  File. 

31  CLOSE  ( UNI T= 1 1 ) 

9999  RETURN 
END 


Open  the  input  data  file 


& 


u 


b 


c 


OPEN  ( UM 1  T=  1 1 ) 
DO  30  K- 1,1000 


C 


Transter  a  record  trom  the  input 
■file  to  the  buffer  SCORE. 


READ  ( 11 , 1 lOl , END=31 )  SCORE 
1101  FORMAT  ( A7’2  ) 

WRITE  (13,1302)  SCORE 
1302  FORMAT  (’  ’,A72) 


C  - If  SCORE  (1:1)  is  a  •>’,  then  SCORE 

C  should  contain  data  and  a  valid 

C  keyword.  So,  see  if  SCURE(2:*ii) 

C  does  contain  a  valid  keyword.  If  it 

C  does,  transfer  the  data  from  SCORE 

C  to  the  input  variable. 

IF  (SCORE <l:l>  -EG.  ’>’>  THEN 

IF  (SCORE (2: A)  .EG.  ’VIE’)  THEN 

READ  (  iC0f>E  <6:  20)  ,  M  1 15)  *  )  OILcflT 
El  SE 

IF  ( SCORE (2:4)  .EQ.  'ROT')  THEN 

READ  (SCORE (o! 20)  ,  ’  ( 1 15)  ’  )  ROTLMT 
ELSE 

IF  (SCORE (2: 4)  .EG.  ’LOL’)  THEN 

READ  (SCORE (S: 20) ,’( 115) ’ )  LOLLMT 
ELSE 

IF  ( SCORE (2:4)  .EU.  ’ HDR ’ )  THEN 
HEADER ( 1 : 30 )  *  SCQRE(6:35> 

ELSE 

IF  ( SCORE (2:4)  .EG.  'DEG')  THEN 

READ  (SCORE (6! 20) , ’ (E15. ?) ’ )  DEuUIL 
ELSE 


L 


C 


10 


- It  :_,CGI<E  (2(2)  is  a  *  Y  ’  ,  then  SCORE 

may  contain  a  Dunham  Coefficient. 
Use  DIGIT  to  find  which  one  it  1  a . 

IF  (SCORE (2: 2)  .EG.  ’  Y  ’  )  THEN 
DO  10  1=0,9 

IF  (SCORE (3: 3)  .EG.  DIGIT ( 1 ) )  THEN 
II  =  I 
END  IF 
CON  TINUE 
DO  20  1=0,9 

IF  (SCGRE(4:4)  .EG.  DIGIT ( I ) )  THEN 
JJ  =  I 
END  I F 
CONTINUE 
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Program:  SUBROUTINE  RDRDUN 

Verbion:  8*1. 0t>'.  13 

Author:  Pau  J  H.  Ostilu-k 

Air  Force  Institute  of  Technology 
Wr i gh t -Patter  son  Air  Force  Base,  OH 


Do  sc  r j  p  t : on : 


This  routine  opens  an  input  1 l 1 e  (unit  11)  and  reads 
all  records  within.  Each  record  read  is  written  to 
the  output  listing  (unit  13).  Data  records  are 
Marked  by  a  ’>’  in  column  1.  These  records  contain 
data  referenced  by  a  single  key  word.  Ail  other 
records  are  considered  comments.  This  routine  uses 
an  internal  read,  eg  READ  (SCORE (SI  19) ,  ’  (EiS. 7)  ’ >  X 
i  eads  rt'OH  columns  5  to  19  uf  the  character 
juriaLle  SCORE  using  the  edit  descriptor  Elb.7  into 
the  real  Variable  X. 


SUBROUTINE  I:  DR  DUN 

CHARACTER* 1  DIGIT (0:9) 
CHAPACTER*30  hEADER 
CHARACTER *72  SCORE 


INTEGER 


1,  II,  J,  JJ,  LVLLMT,  ROTLMT,  U1BLMT 


PEAL  Y(0:v,0:9),  DEOUIL 

70MM0N  /  H1)R  /  HEADER 

COMMON  /DAT/  Y,  DEGUIL,  U I BLMT ,  ROTLMT,  LVLLMT 
DATA  DIGIT  /  ’  0  '  ,  ’  1  ’  ,  ’  2  ’  ,  ’  3  ’  ,  ’  **  ’  ,  ’  5  ’  ,  ’  6  ’  ,  ’  7  '  ,  ’  8  ’  ,  ’  9  ’  / 


WRITE  (13,1301) 

1301  FORMAT  ( / ,  ’ 0 ’  , T  7  1  ,  'Input  Data',/,’  ’,T31,’ 


Initialize  the  input  variables. 


VIBLMT  -  0 
ROTLMT  =■  0 
LVLLMT  =*  0 


DO  S  1-0,9 
DU  4  J-0,9 
Y  (  I  ,  J  )  =  0 
CONTINUE 
CONTINUE 


- Look  tor  (jol  t  Bb. 

Du  20  1-1,26 

IF  iwCUkt  (  3:  -4)  .  LU  .  WUHD26  (  1  )  )  1HLN 

DU  320  J-l , 2's 

it  (  (  I  -  1  )  .  L  U  ■  STrtTk.  (  ILuL  ,  j  )  )  f  FiEN 

CURLUl  -  I  -  1 

- Subr  out  j  rie  Lb'LEFT  returris  the 

absolute  level  number  (LVLNu)  tor  a 
For  a  given  level  relative  to  a 
given  electronic  State. 

CALL  LVLLi-T  iLVLNOi 

Li  -  uV'LUG 
END  IF 
CUNT  ii-iUL 

II-  (LI  .  GX.  SILL)  SILL  =  LI 
Gu  TO  21 
Li  ID  IF 
COl'lT  II'IUL 
oU  i 0  SI 

- Look  For  part  C. 

DO  330  1  =  1 , 10 

lr  (  SCOkL  i  b :  b )  .Du.  Sl  LDL  <  X  >  )  THEN 
C  Ur-.O T  ( 1 ;  i  )  —  SCOkL  ( t»:  b ) 

1  LDL  =  I 

e.JDiF 

CGNTiliUE 

- - - - - - Look  i  or  part  DD. 

DU  SO  i- 1 , 2s 

Xi  <  SCORE  <=»:  ?>  .LU.  NUMDLo(I))  THCW 
Du  3-tC  J  =  l,25 

IF  ((1-1)  .LQ.  STATE  1  I LDL, J ) )  THEN 
CUKLVL  =1-1 

- Subroutine  LVLEl-T  returns  the 

absolute  level  nuinber  (LVLNU)  tor  a 
For  a  given  level  relative  to  a 
given  electronic  state. 

CALL  LVLEF  I  i  LVLNO  ) 

L  2  =  LVcNG 
END  IF 
CONTINUE 

IF'  <L2  . 01  .  SIZE)  S1_L  =  L2 
FOUND  =  .TRUE, 
ou  ro  3i 

EliDi'K 

sGNTINUE 


- - - j-f  u  correct  keyword  has  been  found 

■for  a  transi  t  ion  -from  elate  ’A’, 
level  'BB'  to  state  ’C’,  level  ’bb’ 
Mow  sepat  ate  the  value  for  the 
transition  line  from  its  wei^huri^ 
factoi •  The/  are  separated  by  a  1 9 

IF  (FOUND)  THUN 
FCOL  -  0 
LCGL  =  0 
I- ORI’I  -  ’  (£15.7)  ’ 
to  TO  J-y, 72 

IF  (  SC  L>RE  (  0  I  J  )  .  I’lL. .  '  ‘  .  wMu. 

SCORE (J:J>  . ML .  ’;’)  THEN 

FCOL  ^  J 
CO  TO  41 
CilUlF 
C OH  T i HUE 

IF  (FCOL  .or.  O)  THuN 
Du  tiu  J-FCuL,72 

IF  (SCOREiJiJ)  .Eu.  ’  ’  .OR. 

scure(j:j)  .  eu.  then 

LCOL  -  J  -  1 
Gu  TO  01 
END  I F 
CONI INUE 

FORM  COM)  -  NUI’lBEo  ( LCOL  -  FCOL  1  )  (1:2) 

IF  (  (LCOL-FCOL)  .LT.  6)  FORM  (6:6)  -  FUKI'I  (414) 

- ~ - Fbuiid  the  value  for  the  transition 

1  1  lie  • 

READ  ( SCORE (FCOL : LCOL ) , FORM)  LXNEIL1.L2) 

FCOL  =  0 

FORM  -  ’  (El  5. 7)  ’ 

DO  60  J- (LCQL+ 1 ) , 72 

IF  ( SCORE ( J : j )  . NE.  ’  ’  .AND. 

SCORE  CfJJ)  .NE.  ’!’)  THEN 
FCOL  =  J 
GO  TO  ci  1 
END  IF 
CONTINUE 

IF  (FCOL  .GT.  O)  THEN 
DO  70  J -F  COL , 72 

IF  <SC0RE!J:J)  .EQ.  ’  ’>  THEN 

LCOL  -  J  -  1 

GO  TO  71 
END  IF 
CONTINUE 

FORM  CD: 4 )  -  NUMD26 (LCOu-FCOL  +  1 )  ( 1 :  2) 

IF  ( (LCUL-FCOL)  . LT .  to )  FORM (Gib)  =-  F0RM(4:4) 


Found  the  value  for  the  weighting 
factor. 


kEAD  ( iCOkli  i  FCOU :  LCOL  )  ,  FORM )  WEI GH'I  ( L 1 ,  L2 ) 
ELSE 


Did  not  find  a  value  for  the 
weighting  factor ,  so  it  defaults 
to  1.0 


WEIGHT <L1 , L2>  =  1.0 

END  IF 
END  IF 
END  IF 
END  IF 
END  IF 
EMU  IF 

end;!-' 

ECu  il  I  NUE 


Close  the  input  file. 


CLOGE  l UNIT  —  11) 


ko'l  UKN 


v  - 


c 

C  1-f  :  SUBROUTINE  LVLEFT 

C 

C  Version:  bA.Ub.i3 

C 

C  Author:  Paul  H.  UiiaisK 

C 

C  Air  Farce  institute  or  Techno  logy 

C  vh  lylit-Patterson  Air  Farce  Beee,  Gri 

C. 

C  j«  -iL-  !'■  i  p  t ;  on  :  Ttiiti  routine  returns*  art  absolute  efier  yy  level  riumoer 

C  used  by  LFI  I  "from  the  erict  yy  level  number  CkULVL 

C  that  is  relative  to  tne  electronic  state  CURST. 

L 


subroutine  lvleft  <lvlngi 

CHAkACT  ER*1  CUfoi  ,  STLBL(IO) 

111  TEGER  STATE  (  10, 25)  ,  CURLVL,  X,  J,  LVLNG,  NlNsTUO), 
t  OFFSET  ,  SETCHT 

COMMON  /LABELS/  STLfcL ,  CURST 

COMMON  /  LEVELS  /  N1N2T,  STATE,  CURLVL,  SETCNT 
OFFSET  =  0 

L  - - Matcn  the  current  state  with  a.  state 

C  in  STLBL . 

DO  20  1-1,10 

IF  (brUiL(I)  .EO.  CURST)  THEN 


C  - Now  match  the  current  level  with  a 

C  level  in  STATE ( state , 1 eve  1 ) . 


DO  10  J- 1 , MINST  I i  1 

IF  ;SfATE(i,J)  . EO .  CUkLVL)  THEN 


C  - EVLNQ  is  the  absolute  level  number. 

C  The  level  number  in  CURLVL  is 

C  relative  to  the  state  in  CURST. 


LVLNO  ^  OFFSET  *  J 
GO  TO  9999 
END  IF 

10  CONTINUE 

ELSE 

OFFSET  =  OFFSET  +  NINST ( 1 ) 
LI  ID  IF 

-u  CONTINUE 

9990  RLIUKN 
END 
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INTEGER  ISTP,  IPRINT,  JR,  JG,  JA,  JJ,  EVMCNT,  PARMNO,  NELMT , 
■»  hOEVaL,  JObl'l,  N,  PR,  PP,  EVMLVL  ( 20 )  ,  NGPNTS,  LEVEL, 

*  KM,  RW 
1 14 V  EGER >.  o  I U ,  IX 

REAL  LVM(2b>,  STPGZE,  FARM (10),  PARHLG  (  IO)  ,  PARMHI(IO),  A  (70), 

*  LEG IN ,  EVALS202),  LVEC ( 202 , 202 > ,  H(205G3),  S(202,4), 

<  EVMW125),  RE  1 1  b  (  2t> )  ,  CONST  (10),  NODE  (101), 

+  F  3 1 VAL ( 202 )  ,  HEAR,  MU,  MIN,  L(202,4) 

COMMON  /ENERGY/  EVM,  EVMW,  RESID,  EVI'ILVL,  EVMCNT,  STPS2E, 

+  NELMT,  bEGIN,  HEAR,  MU 

COMMON  /FARMS/  PARI-I,  PARMLG,  PARMH1,  CONST,  PARMNO 
COMMON  /EIGENS/  EVAL,  EVEC,  NuEVAL,  JOBN 
Col-IMUN  /MATRIX/  H,  S,  L,  N 

COMMON  /WAVFUN/  NODE,  PSIVAL,  NORNTS,  LEVEL 
COMMON  / CMRLEL /  STLKL 


e:;ter mal  pun 


-Vtii  the  uutput  listing  tile,  print 
the  header,  end  &tai t  run  stats. 


C  it L  1 1»_  n Lu k 


.  - - - Head  the  input  tiles  tor  control 

parameters,  input  data,  arid  random 
number  gener  ator  bcedb. 

CnLL  READCk  (IU,  l'..,  IS  TP,  IPr,  li'IT ,  Jk,  JG,  JA,  JJ,  PR,  PP ,  Rl<1,  RW 
IF  \  KM  .  NE .  0.'  THEN 


Pirid  the  aet  ot  potential  energy 
parameters  that  gives  the  best 
least  squares  tit  ot  measured 
and  calculated  energy  values. 


3  OBI  l  ~  0 


CALL  WI  MUM  ( F ARMNO , HARM , A, FUN,  iU,  IX,  IlTP,  I PR 1 NT , J R , j G , J A , J J > 


- - Put  the  current  random  number 

generator  seeds  into  a  tile  tor  use 
on  the  next  execution  ot  this 
program . 


Call  putrwi)  (iu,  ix) 
till)  ip 

;l:  (RW  .  NE .  0)  THEN 


Find  the  normalized  nave  tunctioris 
that  correspond  with  the  calculateu 
energy  levels. 


JOHN  -  1 
CALL  WAVE 


HIM  -  FUI-HPAKM) 


LNDiP 


Create  the  output  listing  and  tiles 
with  residual  and  potential  plot 
data,  and  the  wave  functions. 


CALL  OUTPUT  <Pk.  PP) 


ClOiif  the  Output  ♦  1  1  t  OHIO 
ututjstico. 


CALL  i Km  I LP 


9VV9  CPU 
4>AI)Q  t  HLADLK 
iriDL,  KLADEK 
TAblf,  l-U  i  KND 
i-ADD,  FUN 
Wibb , POTENT 

Fmi)J3  ,  L  I  OLl'l 
j>m£D , C  LTL 
iAjj, F0LD2 
•P  A  L  D ,  FOLDX 
*mPD  ,  UAl'C 
+  mL-L  .  URFGLD 
i  .-il  i) ,  M  L1  r\  rl  A  L 
It  l.-V,  OUTPUT 
ltd,  C'ETF’o  r 
TnL D ,  PkilTLK 
i  .-.Li, ,  t  LTRLi 
I'A.L’P  ,  PL  I  P  O  i 
-i  .-.00  ,  PU  i  vJAU 
i  iuiU,  I  P.A  I  LP. 


6 
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C  Program:  SUEkOUTiNE  HEADER 

C 

C  Vcriitin:  o-( .  1  1  .  -  u 

C 

C  rtuttur:  Puuil  H.  (JL.Uiek 

Air  Farce  Inutitute  ol  TecLr.o  i  u‘J ./ 

C  Ur  i  jh  t -Pat  t  c*r  i-on  Air  Force  Lcbc,  GH 

C 

C  Geeci  iption:  rtLYiEER  ojjei..  trie  output  iibting  file  (  logical  unitb 

1  i3  and  a)  af.J  ituft  j  CPU  and  wall  tune  uae 

'  btatiiitlCa. 

C 

SUEKOU  TINE  HEADER 

CHAkACTEP.*8  VERSN 
CHAkaCTER*13  PCM 

INTEGER*3  I T)  A  T  E  (  3  )  ,  ITII'ltO),  IUSERW) 


"t'P  GM  *  ’  84 .11. 30  ’ 

•CM  -  ’ 0EP/o4D-o/ 1 . 3 ' 

C  - li.it  late  the  run  statistic®. 


Cr\LL  ETIME 
CALL  STINE 

C  -  - - - - - Get  the  current  duU,  time,  and 

C  ucer  tidi.ic  ror  output  on  the  header 

-hLL  Ln i  L  .  .  m i l i 

C  r  .LL  TiliciZr.liE) 

0  ALL  ■  •  .  CP  i  iu  (  1  US  Li-: ) 

- - - - Open  the  output  listing  Pile  ana 

•1  write  out  the  header 

OPEN  ( UN  I T - 1 G ) 

OPEN  (Uttir-d) 

MITE  ill,  1301)  lUuER,  VERSN,  IDATE,  FCH,  lilHu 
130  1  FORMAT  i’i  liter  :  ’  ,  4AG  ,  Til ,  ’  A  i  r  Force  Institute  ot  I ec hno 1 oyy ' , 

•  T  1  1  0,  ’  Ver  ia  1  on  :  ’  ,  lio, 

r  /,'  Date!  ’  ,  G  A  G  ,  T 1 1  A  ,  ’  F’  C  M  l  *  ,  A  1 G  , 

♦  Til..--:  ’  ,  3  AG ,  T  64  ,  ’  D I A  T  ON  ’,/,/) 

nr,,.-,  (.  L'-UPil 

END 

F-4 


P  •  ugr  am 


K  Lm  L:  L  K 


Vcf  s  l  o  n  :  o  4  .  11.  jO 
Author:  Paul  H.  Ostdit'k 

Air  Force  Institute  o)  T ec hno 1 oyy 

Wr  i  <jh  t  -Ps  t  ter  son  Air  Force  Ease,  OH 

Dtr  s  c  r  i  p  t  i  on :  READER  opens  3  input  tiles  (unit  11,  12,  13)  and  reads 

el l  records  within.  Each  record  read  is  written  to 
tits*  output  listing  (unit  13).  Data  records  are 
marked  by  a  in  column  1.  These  records  contain 

data  r  of  er  enced  by  a  single  key  word.  Ail  other 
records  are  considered  cowmen  cs .  This  routine  uses 
an  internal  read,  ey  :<LnLi  ( SCORE  <  ‘a :  1  / )  ,  '  ( E  1  b .  7 )  '  )  X 
l  eads  i  r  aiti  columns  t>  to  19  of  the  char  acter 
. ar  ruble  sC&KE  usiny  the  edit  descriptor  els. 7  into 
t  )'t  a  r'esl  V  al‘  1  aO  1  e  A  • 


eGBROU  T  X  HE  KLnUtk  HU,  IK,  IETF,  1  PRINT  ,  JR  ,  J  0  ,  Jn.JJ  ,  PR  ,  PP  ,  KM ,  KW  ) 


3 ;  rAKCAC TER*  1 
ChAr<ACTER*2 
s.-iARACT£R*9 


DIGIT ( 0 : V ) 

LOCI'.  IV),  PLUCK, 
FORM 


CMARAC  fliKKs’S  ECUr.E,  3  i  LltL 


NUMBER  i  0  1  30 ) 


1  iJ  IcGEK  LuOR!4/,  t-VMCNT,  I,  IRKIN'),  iolP,  J,  jh,  PR,  PP, 

JG,  JJ,  JR,  IJELH  i  ,  PARl'itiu,  F COUNT,  FCGL,  ECOE, 

E VI'ILOL  l.  jt  ,  Ail,  KW 

1  in  I  LutKIO  It),  I/. 

LdOl CmL  L  UK l LI 

REAL  EVMdl),  .-mkI'I  (  10)  ,  PARMLO(lu),  PARMhl(iO),  STrSZE,  BEGIN, 
ENDS,  TEMP,  REo iD ( 3j ) ,  EVMW(25),  CONST (10),  HBAR,  MU 

i'oMMON  /ENERGY/  EVM,  EVMW,  RESID,  EVMLVL,  EVMCNT,  SI’PGZE, 

NELM  T ,  BEG  IN,  HEAR,  MU 

COMMON  /FARM a/  FARM,  i'ARMLU,  PARMHI  ,  CONST,  PARMNO 
COMMON  / CHPLDL /  3TLBL 


Data 

DIGIT  /•O’,’ 

i 

t  » 

9 

»  »  •  i 

»  w 

*  4  *  »  *=■  »  * 

9  ^  9  9 

o’  ,  '7 

c: 

nC 

DAT  A 

NUMBER  /’GO’ 

$ 

’  01  ’ 

,  ’02’ 

,  ’03’ ,  ’03’ 

,  ’03’ 

, ’Oq’ , ’07’ 

> 

* 

00  ’  , 

’09’  , 

’  10 • ,  *  1  I  ’ 

f 

'  12  ’ 

,  ’  13' 

'  It,’ 

,  ’  16’ 

,’17’,’IG’ 

9 

19’  , 

’20’  , 

'31  1 

9 

1  —  1 

9  - 

,  ’  14 ’  ,  ’ 2o ’ 

,  '  2o  ’ 

, ’27’ , ’20' 

9 

2y  ’  , 

’  30  ’  / 

DA  I  A 

LOCK  / ’ NL ’ , ’ 

I 

> 

II-  '  ,  ■ 

JR’  ,  ’ JO"  ,  ’ 

JA’  ,  ’ 

J J ’ , ’ PR’ , ’ 

hP 

■  / 

—  Irijit.iij.iizc.'  ca  o  in  e  uai'laulea 


b  0  10  i  -  i  , 

I  ui 0 i ■  i  i  >  —  <-• 

i  0  0  L'l  J  I  1  (  JOE 

I  nAMItl.1  —  U 

C  - - - - Open  the  Control  Parameter  Input 

C  Pile  --  logical  urilt  li. 


GPU)  lUHir-ili 
v:  i  lie.  i  i  j  ,  i  j  1  / 

1  .  1  TOP:  il  AT  (' O’  ,T24,’ Control  Pat  uirie  ter  input  Pile’,/, 


DO  AO  1-1,1 000 


C  - - - Vi  oMblc-r  a  record  from  the  input 

C  file  to  t  lie  buffer  SCORE. 

READ  ( II , I  101 , LNb  =  4l )  SCORE 
lie  I  FORMAT  (AS2) 

WRITE  (13, 1-02)  SCORE 
1-02  F  ORl'tA  C  (’  ',a:'21 

C  - - - - It  SCOrcE  (  X  :  i  )  lb  a  ’  /  ’  ,  then  SCORE 

C  bheuid  contain  data  and  a  valid 

C  keyword.  So,  bee  it  SCGRb.(2:s) 

doeb  contain  a  valid  keyword.  If  it 
0  dues,  transfer  the  data  from  SCORE 

C  to  the  input  variable. 


IF  (SCORE  till)  .  Eu.  ’>’)  THEN 
Lr  ( SCORE (2:2)  ,£u.  ’hB’)  THEM 

READ  < SCORE (b: 1 V) , ’ (Eib. 7) ’ >  HbnR 
ELSE 

IF-  <  SCORE  (l:  2)  .  El«.  ’HU’)  THEN 

r  :LmU  (  j2uKE  (  j(  i  V  >  ,  ’  IE1  j.  7)  ’  )  MU 

EL.  sE 

IF  ( -uCOkL  (  2 1  3 )  .EQ.  ’  bO ’  )  THEN 

READ  (SCORE  4b: 1V> ,  ’  (Elb. 7)  ’  >  EEC  I N 
ELSE 

IF  ( SC  UKl ( 2  I  2  J  .Eu.  ’EN’)  THEN 

RE mD  ioCGKE 4b; IV) , *  (CIS. 7) ’ )  ENDS 
ELSE 

IF  luuoiE(-;JI  .EU.  ’KM’)  THEN 
hu,-,D  (SCORE  (b:  19)  ,  ’  (  I  ib)  ’  )  KiT 

LLL.il 

II  (SCORE (2: 2)  . EO .  ’KW’)  THEN 

READ  ( SCURL (b;  1‘7),  ’(lib)’)  RW 
EL  .  L 
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j  I  >  jLwt-.u  i  2  1  2  ) 

DO  20  J  -  i  ,  1  0 


’  t-  )  i  I1LN 


*  Lu. 


ivLCui  i  i :  l )  -  ’ » ■  ■ 

IXuCI-  (2:2)  -  NUHbLk  ( J  )  <2:  2) 

it  j)  .tu.  KLuCKi  I  rli_w 

KlAU  (  Sv-OkE  ( E>  1  i  9  )  ,  ’  l  c.  1  3  .  /  )  ’  )  F  Aim'I  i  J  ) 
ri-ikl'lll  u  —  FtiRNHu  i 

uu  ro  -t  o 

END  IF 
CONI iNUC 
tLL-t 

It  (s^GKE  (2:  2)  .  Lu.  *  U’  )  THEN 

JO  CL  J-l , 10 

(.LOCK  l  1  :  1  )  -  'O' 

FLOCK  (2:2)  -  NUMoCR ( J ) (2: 2) 

IF  (SCORE  (2:  3 )  .  C&.  RLOCK)  THEN 

IstnO  (SCORE  (S:  19)  ,  ’  iE13.  7)  ’  )  FARNriilJ) 

00  TO  AO 
END  IF 
C  OUT  i  NUC 
ELSE 

IF  (SCORE  (2:  2)  .EC;.  ’L’)  THEN 

JO  22  J-l , 10 

KEOCK(i:l>  -  ’E’ 

KLUCK  (2:2)  -  NUilbt_R  ( J  )  <2:  2) 

it;  (SCORE  (2:  S>  .  Cu.  RLOCK)  I  HEN 

READ  (SCORE  (Si  l'>)  ,  ’  (E15.  V)  ‘  )  PARI'lLO(J) 

0 U  TO  AO 
LIIUIF 
C JN , 1UUL 
ELSE 

IF  (SCORE (21 2)  . CO.  ’C’>  THEN 

DO  23  J-l , 10 

FLOCK ( 1 : 1 )  =  ’ C ’ 

(•.LOCK  (2)  2)  -  NUITLLR  ( J  )  (2:2) 

IF  i SCORE (2:3)  .  EU.  RLOCK)  THEN 

READ  (SC0KE(b:19) , ’ (E1S.7) ’ )  CONaT(J) 
GO  10  So 
Li  Wit 

cunt inul 

ELSE 

DO  30  K- 1 , 9 

IF  (SCORE (2: 2)  .  EQ.  LOCK(K))  THEN 

Ki_nD  (SCORE  (S:  lV)  ,’(  113)  ’  )  DOOR  (K) 

EUs  I F 
Ct-A-iT  1  UUE 
END  IF 
END  IF 
Li  4 is  I  F 
END  11 
E  N  J  i  t" 

END  1 1 
EN2  i  I 


1.1:0  1 1 


C  Program:  SUBROUTINE  UNFOLD 

C 

C  Ver  s.  l  or :  03.07.84 

C 

C  Author:  Paul  li.  Ostdicl. 

C 

C  Air  Force  Institute  of  Technology 

C  Wr  1  gh  t -Pat:  ter  son  Air  Force  liabe,  OH 

C 

C  Description:  UNFOLD  recovers  the  eigenvectors  v  of  Hv-eSv=0 

C  from  the  eigenvectors  y  of  X/-ey=0  where  e  are 

C  eigenvalues. 

C 

c - 

SUBROUTINE  UNFOLD 
INTEGER  I,  3,  K,  N,  NOEVAL 

REAL  E UAL  t  202 )  ,  EUEC ( 202 , 202 )  ,  HC20S03), 

+  8(202,4),  SUN,  L ( 202 , 4 ) 

COMMON  /EIGENS/  EVAL,  EOEC  ,  N0EUAL,  JOBN 
COMMON  /MATRIX/  H,  S,  L,  N 

LO  30  I -N , 1,-1 

DO  20  i - 1 , NOEVAL 
SUM  -  0 
DO  10  K-1,3 

IF  (  (  I  rf  )  . LL.  N)  THEN 

SUM  -  SUM  ♦  <L(  (  I  ♦  K  )  ,  (4-K)  )  *  EVEC ( ( I +K ) , J )  ) 

END  IF 

10  CONTINUE 

EVECtl.J)  -  ( E  0  £  C ( I , J )  -  SUM)  /  L(I,4) 

20  C  1NTINUE 

0  CONTINUE 

Pi  I  URN 
lNU 


('r  L.^11  vtl.l 


SUtKOLiTINL  FOLi):< 


t  N I EGER  I,  J,  K,  KK,  LL,  LLL,  N 

REAL  H(20-dG3),  L<202,4),  3(202,4),  SUM 

COMMON  /MAT  KIN/  H,  3,  L,  N 


ito  :.u  j  =  i,n 

b‘J  2  0  1  -  J  ,  N 
SUM  -  0 
1)0  10  K- 1,3 
LL  =  J  -  K 
LLL  =  4  -  K 
IP  (LL  .  GT.  0)  THEM 

kk  =  (I  *  <  I  -  1  )  /  2)  +  LL 

SUM  -  SUM  +  <H< KK)  k  l(J,LLL)) 

END  IF 

i  '  CUM  I  I HUE 

K  =U*(l-i)/2)+J 
MU)  -  (H(K)  -  SUM)  /  L  ( J ,  4 ) 

20  CONTINUE 

30  CONTINUE 


9o  (•  L  TURN 

r  md 
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c 


C  Hr  LKH  u. .t:  'jULKL'Ij  !  -NE  FOLDS 

C  Vcr  l  -if;  :  04  .11. 30 

C 

C  rtathur:  L-aul  H.  0-.t<iieL 

C 

C  Air  Force  Institute  at  i  cL  h  f  ill,  i  Oy  / 

C  w'r  i  gh  t -Pat  tfi'aon  Air  Furci  Laae,  OH 

C 

C  Dc&c  r  i  p  t  i  ori :  T 1. 1  l.  icutimr  t  inas  a  matrix  Z  sucli  that  H  -  Li. 

C  Nutr  ix  Z  is  stored  in  pUce  oi  H  l  ri  the  array  H. 


SUEROUT I  HE  FOLDS 

IHi'LOLR  CNT,  1,  II,  J,  K,  KK,  LL,  Id 

RdhL  H(20303>,  H3E  ( 202 , 3 )  ,  L(2G2,4>,  S(202,4),  SUM 
COMMON  /MATklX/  H,  3,  L,  Id 


DO  30  I  -  1  ,  W 

1 F  (  (1*3)  .  LT.  N>  1  HON 

CNT  *  I  ♦  3 
EL'uE 

CNT'  =  N 
END!  F 

DO  20  J - 1 , CNT 
-  0 


DO  IO  K-l, 

■r 

LL  -  4  - 

K 

I  I  -  I 

K 

IF  (11  . 

OT 

.  0  > 

THEN 

IF  (II 

. 

GE  . 

J )  THEN 

KK 

= 

(  I  I 

*  (11-1)  /  2)  +  J 

SUN 

-= 

SUM 

+  ( L ( I , LL )  *  H ( KK ) ) 

ELSE 

IF  ( 

v  J 

-  1  I  i 

.LL.  _•)  THEN 

CUM  =  SON  +  ( L ( I , LL )  *  HSU ( I  I ,  (  J- I  I  )  )  ) 

LliDIP 
END  IF 
End  i 

i  J  0  Or.  I  I  I  IL  L. 

. L  (I  . 00 .  j)  f  HEN 

K  -  ti  4  t  X  —  1  )  /  2)  +  J 

H ( K >  -  ( H  <  K )  -  SUN )  /  L (  I  , A ) 

ELsE 

K  ~(J*(J-l)/2)+l 

H3B(  I  ,  (J-I)  >  =  (H(K)  -  SUIT)  /  HI, 4) 

END  IF 

LO  CONTINUE 

30  CONTINUE 

'->999  RETURN 
END 
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SUBROUTINE  GETL 


Vefcieirtl  U4.  1  1  .  _>0 

nut)iOi  :  Paul  H.  Getciiel. 

Mir  Force  Institute  oi  Technology 
Ur  i<jtit-i  attersan  Air  Force  Baee,  GH 

0  ••  ■ .  i .  t  lpt  ior.:  This  r  out  me  dec  o=,e=>  the  matrix  S  into  a  matrix  l. 

•juch  that  ‘S  —  LL(transpOse)  . 


SUBROUTINE  Ot  TL 

INTLGLK  I,  II,  ILI’IT,  J  ,  LL,  N 

PlmL  H ( 2GSG2 ) ,  3(202,4),  SUN,  L l 202 , 4 ) 

COMMON  /WmTRI/'/  PI,  3,  L,  N 

I L  NT  ^  4 

DO  20  I  -  1  ,  h 

UO  10  X  —  J  ,  I  LIT  r 
SUIT  -  0 

. F  (l  .EG.  J)  THEN 

SUM  a  L<I,l/**2  ♦  L(I,2)**2  +  L (I , 3 ) **2 
L ( I , 4 )  =  SORT  IS (1,4)  -  SUM) 

ELSE 

IL  -  4  -  I  ♦  J 
IP  (LL  .Lu.  2)  THEN 

SUN  =  L  (  I  ,  1  )  *  L(J,3> 

ELSE 

IP  (LL  .EG.  Z)  THEN 

SUN  =  L(I,1)*L(J,2)  *  L(I,2)KL(J,3) 

END  IP 
ENDIP 

L ( 1 , LL )  =  (3(1, LL )  -  SUM)  /  L(J,4) 

PI NL  I P 
CUNI i HUE 

i I  -  MUD ( J , 2 ) 

IP  (II  .EG.  0)  THEN 
ILMT  -  J  »  4 
EL'  L 

ILMT  -  J  ♦  3 
END  If 
CONTINUE 


RETURN 


n  n 


CALL  FOLDX 


L ■  M> 
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n  n  n  n  n  r>  n  n  n  o  non  o  n  r>  n  n  n  n  o  n  n  n  n  n  n  n  n 


Pr  Ci'jr  jl:ii  :  SUBROUT  I  HU  EIGEN 

Vc-r  -  i  u-ri :  UA  .11. 30 

Author:  Paul  H.  GstdieL 

Air  Force  Institute  of  Technology 
Wright-Pattcrson  Air  Force  Base,  OH 


Descr  i  p  t  i  on  : 


This  routine  -finds  the  eigenvalues  and  optionally  the 
eigenvec  tore  <JQBN=i)  o-f  the  generalized  eigenvalue 
problem  Hv-eSv-0.  This  is  done  using  routines  GETL, 
F0LD2 ,  FOLD;:,  UNFOLD,  and  IMSL  routines  EHOUSS, 
EQRT1S,  EGRT2S,  EHOBKG. 


SUBROUTINE  EIGEN 

INTEGER  I,  IER,  ISW,  1Z,  JOBN,  N,  NQEVAL,  M,  Ml 
LOui CAL  FIRST 

REAL  EOAL ( 202 ) ,  EVEC ( 202 , 202 > ,  H120503), 
r  S(202,4),  L(202,4),  SUM,  E(202),  E21202) 

COMMON  /El GENS /  EUAL ,  EVEC ,  NOEVAL,  JOBN 
COMMON  /MATRIX/  H,  S,  L,  N 

OhTh  FIRST  / .TRUE. / 


-Use  the  Cholesky  decomposition  to 
get  a  matrix  L  such 
that:  S  =  LL ( transpose ) 

This  only  has  to  be  done  once  since 
S  doesn’t  change. 


IF  (FIRST)  THEN 
CALL  GETL 
FIRST  =  .FALSE. 
END  IF 


-Get  Z  such  that:  H  -  L2 


CALL  FOLDZ 


-Get  X  such  that:  2  =  XL ( transpose) 

the  problem:  Hv  -  esv  =  0 

now  Becomes:  Xy  -  ey  —  0 

where:  e  —  eigenvalues 

v  -  eigenvectors 
y  -  eigenvectors  such  that: 
y  -  L ( transpose ) v 


1)0  1X0  I-l,ELhlT 

II-'  (EVI'KI)  .EG.  0.0)  THEN 
RESID(I)  =  O 
ELSE 

kESIDil)  -  EVM(X)  -  EOAL(l) 

FUN  =  FUN  +  (EVNW(I)  *  RES  I D  (  X  )  *-k2  )  *  0.5 

END  IF 

O  CONTINUE 

IF  (FUN  .GT.  IE  1 0 )  FUN  =  1E10 

■99  RETURN 
END 


JJLESS  -  0 
DO  30  J- l , 1 0 


Use  the  potential  values  ana  matrix 
b  CL'iia  to  build  the  sub  ~  Ilia  1 1  icies. 


vi  <j  > 

V2  ( J  ) 
V3  (  J  ) 
V  4  (  J  ) 


P0  *  Villi  J) 

!3*P0  «■  DP0*STPG2E>  K  V2H ( J > 
PI  *  V2H ( J ) 

(3*P1  -  1)P1*STPC2E)  *  V4H  <  J  ) 


KK  -  (II  *  I  I  I  - 1  )  /  2)  ♦  JJ 

LL  -  4  -  II  +  JJ 


Load  the  ^uu-iiiali u  values  into  H 
arid  0. 


H(KK)  =  H ( KK )  +  I ( j )  ♦  VI  (J>  +  V2(J)  +  V3(J)  ♦  V4(J) 

IP  (FIRST)  THEN 

Sill, LL )  -  S(I1,LL)  +  SH ( J ) 

Li. DIF 

i.  (II  .  L  tt .  J  3  )  THtlN 
II  ^  II  +  1 

JJ  =■'  JJ  -  J  J  LESS 

JJLESS  -  JJLESS  +  1 
ELSE 

J  J  -  J  J  +  1 
END  IF 

i:  Si'!  r  i  nue 

II  =  II  -  2 
JJ  --  JJ  +  2 
CC NT  I NUE 


FIRST  —  .FALSE. 


Subroutine  EIGEN  solves  the  general 
eigenvalue  problem; 


HV  - 

eSv  *  0 

where!  H 

1  ^ 

tlie 

H  (energy )  matrix 

v 

2  & 

t  he 

e  2  gen vec  tor 

e 

l  * 

the 

eigenvalue 

c? 

i  & 

the 

normal iiing  matrix 

lAlL  eigen 

FUN  =  0.0 
LLl'ir  -  EVMCNT 

IP  (NGEVAL  -LT.  EVt'ICNT )  ELMT  =  NOEVnL 

- - -Compute  the  residuals  and  the  value 

o+  FUN. 


P-14 


V  til  (  3) 

2uG0 

* 

S  T  P  s  2  L  *  *  2 

/ 

3c,2ts0-_. 

VI r:  <  A  » 

SOAUU 

* 

^  TP GEE 

/ 

oouO 

V  I If  it.) 

iJVvjlJ 

* 

Si  PELL  *.  A  2 

/ 

J>  u  If  b  C  0  0 

V  A  i  i  (3 

-  lG->c~0 

* 

GTP3LE 

/ 

362CG00 

v-.i  ■:  (’ ) 

-  -  1  1  ;2o 

k 

STPGLL  a  X2 

/ 

Z&lz  dbOO 

V-li :  ;  •  - 

-  -20 JO 

k 

~i  1  f  ^  O 

/ 

o  b  iJUU 

v  li :  i  ) 

=-  -2. .0-10 

* 

S1P  .  2E  a-  k 2 

/ 

3612  bbOO 

VAH ; 10) 

—  i  u  L' 

k 

3TFS2L**3 

/ 

3628606 

- - - - - - N  1  s  the  number  of  members  on  the 

mean  diagonals  of  H  and  S. 

N  -  (2  X  IIELMV )  *  2 

l  NO  IF 

- - -  -  -  - Clear  out  the  0rra/b  that  H  arid  S 

are  stored  in.  Do  this  only  the 
first  t  i  me  -for  S. 

DO  20  J -  1  ,  M 
DO  10  1-3, N 

f'K  —  (I  *•  C I  —  1  )  /  2)  ♦  J 

LL  =  -4  -  1  +  J 

HCKK)  -  0 

IF  (FIRST  .AND.  LL  .  GT.  0)  THEM 
3 ( I , LL )  -  0 

ENDIF 
CONTINUE 

CONTINUE 

- - - - - Duild  the  sub-watr i c l es  that  are 

pieced  into  matneiea  H  and  S 
from  the  seeds.  This  is  done, 
in  c*r.  iterative  wanner  ,  once  for 
each  step  along  the  grid. 

f STLP  =  .  M:JE. 

Du  AO  .  - 1  ,  NELI’ll 

f:o  -  i ec ii-i  ♦  (stpgie  *  (i-m 
R 1  -=  ELGIN  +  (GTPGZE  *  I) 

- - Subroutine  POTENT  returns  the 

value  of,  arid  slope  of  the  model 
potential  function  at  each  end  of 
the  current  grid  element  (R0,R1). 

CALL  POTENT  (R0,  Rl,  P0,  DPO,  PI,  DPI) 

IT*  (FSTLP)  THEN 
II  ^  1 

J  J  =1 

FSTLP  -  .FALSE. 

ENDIF 


Program:  FUNCTION  FUN 


C  Version:  64.11.30 

r 

C  Author:  Paul  H.  Gi>tdiek 

C 

C  Air  Foret  Institute  o+  Technology 

C  Ur  1 ght -Pet ter son  Air  Force  Base,  OH 

0 

C  Description:  Tins  routine  solves  the  Schrodinger  nave  equation 

C  using  o  -et  of  potential  energy  parameters  to  build 

G  the  Mave  equation  in  matrix  form.  Then  EIGEN  is 

C  used  tu  soive  the  eigenvalue  problem  for  the  eigeri- 

C  ...lues,.  These  eigenvalues  are  compared  in  a  weight- 

G  td  least  squares  sense  with  the  observed  energy 

C  levels.  The  value  of  rUM  returned  is  the  sum  of  the 

C  Weighted  residuals  squared. 

C 
C 


REAL  FUNCTION  FUN  ( PARM ) 

INTEGER  EVMCNT,  I,  IER,  II,  J  ,  JJ,  JJLESS,  JOHN,  KK ,  LL,  N , 
r  NGEVAL,  MELNT,  ELMT,  EVMLVL (25) 

LOGICAL  FIRST,  FSTLP 

REAL  E  VM ( 25 )  ,  T(10),  VIH(IO),  V2H(10),  V3H(10),  V4H(10),  VI (101, 

+  V  2  (  1 0 )  ,  V3(10),  V4U0),  STPSZE,  R0,  Rl,  H(20503),  S(202,4), 

i-  SHI10),  P0,  DPO,  PI,  DPI,  FARM  (10),  EVAL(202>,  EVEC  ( 202 ,  202 1 

+  DELTA,  RES  I D  ( 25 )  ,  BEGIN,  EVMW(25),  HBAR ,  MU,  PKE ,  L(202,4> 

COMMON  /ENERGY/  E VM ,  EVHW,  RESID,  EVMLVL,  EVMCNT,  STPSZE, 

^  NELMY,  BEGIN,  HBAR,  MU 

COMMON  /El GENS/  EVAL ,  EVEC,  NOEVAL,  JO BN 

COMMON  /MATRIX/  H,  S,  L,  N 

DATA  FIRST  /.TRUE./ 


C 

c 

c 


Load  the  arrays  used  to  build  the 
matricies  H,  S,  and  NORM.  These 
will  not  change,  do  this  once  only. 


(FIRST) 

THEN 

SH  ( 1  ) 

= 

1872 

* 

SH  (2) 

= 

264 

* 

SH  (3) 

= 

48 

* 

SH  ( 4  ) 

= 

648 

* 

EH  (5) 

= 

156 

* 

SH  ( a ) 

= 

1872 

* 

SHI?) 

- 

-lto 

* 

SH  ( S  ) 

- 

~  3  cj 

* 

STPSZE 

/ 

5040 

STPSZE**2 

/ 

5040 

STPSZE**3 

/ 

5040 

STPSZE 

/ 

5040 

STPSZE**2 

/ 

5040 

STPSZE 

/ 

5040 

STPGZE**2 

/ 

5040 

STPSZE**3 

/ 

5040 

END  IF 
LI  ID  1  F 
:CN I 1MUE 


END  IF 
LIlDlF 


L:iJ  1  r 

C  ,'N  I  J  HUE 


-ClDbe  the  input  tile. 


clo^l  ;uhit=12> 


-Open  the  Random  Number  Generator 
Seed  tile  --  logical  unit  1-*. 


OPEN  ( UN I T= 1 4 ) 

WRITE  (13,1304) 

FORMAT  (’ 0’ , T21 ,’ Random  Number  Generator  Seed  File’,/, 
’  ’  ,  T21  ,’ - -  - - - ’) 

DO  120  I  =-1,1000 


-Transfer  a  record  from  the  input 
file  to  the  butter  SCORE. 


READ  (14, 1401, END-121)  SCORE 
FORMAT  ( A72 ) 


WRITE  (13,1302)  SCORE 


-it  SCGRE(i:i)  is  a  tneri  SCORE 

should  contain  data  and  a  valid 
keyword.  So,  see  it  SCGR£<2:3) 
does  contain  a  valid  keyword.  It 
does,  transfer  the  data  from  SCORE 
to  the  input  variable. 


IF  (SCORE (111)  .EQ.  ’>’)  THEN 
1 F  ( SCORE (2:3)  . EU .  ' IU ' )  THEN 

READ  (SCOPE tt:  IV)  ,  ’  \ I  lb)  ’  )  IU 
ELSE 

IF  (..CORE  (2:  3)  . EU.  ’IX’)  THEN 

READ  (SCORE (b: 19) , ’ (lib) * )  IX 
END  IF 
END  IF 
END  J  F; 

C  j.  T  J  i'iUE 


-Close  tlie  inr>  .  file. 


CLOSE  (UNI T- 14 ) 


RETURN 


c 

c 

c 


does  contain  a  valid  keyword 
does,  transfer  the  data  from 
to  the  input  variable. 


at 

31 


>0 

71 


i 

.  1 


If  <SCGRE(l:i)  .EU.  *  >•)  THEN 

IF  (SCORE  (2  : -4)  .  EG.  1  LElL  ’  )  THEN 

STLBL (1:72)  =  SCORE (o: 77 i 

ELSE 

IF  (SCORE  (2:  2)  .EG.  ’V’)  THEN 
DO  210  K-0,25 

IF  (SCORE(3:m  .Eu.  NUMBER  ( i< )  )  THEN 
EVMCNT  -  EVMCNT  t  1 
EVMLVL (EVMCNT)  -  K 
F COL  -  0 
LCOL  =  0 
FORM  =■  ’  (Ei5.  7)  ' 

DO  EG  J -o , 72 

IF  (3CGRE(J:J)  .NE.  ’  *  .AND. 

+  SC0RE(J:J)  .NE.  *;•>  THEN 

FCOL  -  J 
GO  TO  El 
El  ID  IF 
CONI INUE 

IF  (FCOL  .01.  0)  THEN 
DO  70  J-FCOL , 72 

IF  (SCORE(JJJ)  .EU.  ’  *  .OR. 

+  SCORE ( J ; J )  .EQ.  ’ 5 ’ )  THEN 

LCOL  =  J  -  i 
GO  TO  71 
END  IF 
CONTINUE 

FORM (3:4)  =  NUMBER l LCOL -FCOL + 1 > (1:2) 

IF  i  ( LCOL -FCOL  1  .LT.  6)  F0RM(6:6)  =  F0RM(4:*fl 
READ  (SCORE (FCOL: LCOL) , FORM)  EVM (EVMCNT) 

FCOL  =  O 

FORM  =  ME  15. 7)  ’ 

DO  30  J=  (LCOL+l )  , 72 

IF  (3C0RE(J:J)  .NE.  ’  ’  .AND. 

f  SCORE(J:J)  .NE.  ’!’)  THEN 

FCOL  =  J 
GO  TO  SI 
END  I  F 
CONI INUE 

IF  (FCOL  .GT.  0)  THEN 
DO  90  J-FCOL, 72 

IF  (SCORE (J;J>  .EQ.  ’  ’)  THEN 

LCOL  -  J  -  1 
GO  TO  VI 
END  IF 
CONTINUE 

lUkl-li  JI4I  -  NUMBER  (eCUL-FCOL  +  1  )( 1 :  2) 

IF  ((LCuLFCuL)  .LT.  6)  FORM  (  o :  6  )  -  FORM  ( *4 : 
READ  MCOKE (FCOL: LCUL) , FORM)  EVMW(EVMCNi) 
EL'.,L 

EVMW( EVMCNT)  =  1.0 
END  I  F 

F-9 


If  l 

SCORE 


END  »  i 


ENDil 
LlU  i  I 

**i  0  i  IJ.I.  i  '.Ui- 


Clu^e  tilt--  input  file* 


41 


-I  L. 


CLO  -L  (  UN  I T  —  1  1  ) 


NELMT 
1  £  TP 
X  PR I  NT 
JR 
j  G 
J  A 
J  J 
PR 
PP 


DOUR  (  1  ) 
DOOR (2) 
DuUR  (  3 ) 
HOOK l 4 ) 
DOOR  lb) 
DOOR ( a ) 
DOUR  (?) 
DOOR (8) 
DOOR (9) 


(ENDS  - 
Du  4ti  I  —  1  ,  PARHHO 


DLL  III) 


PARl'lHI(i)  =  PARi'lH  1  (  I  ) 
PARI’I(I)  =  (PARMU) 
CONTINUE 


/  HELI'IT 


-  PARi'lLO  <  I  ) 
PAK'MLO  (  I  )  > 


PARMHI ( I ) 


DO  £0  I  -  1 , 2t> 
LVi'l(l)  -  0.0 
EUMLVL(l)  =  0.0 
‘.vu  CONTINUE 


LUMEN  I  -  O 


Open  the  Hfasur  cd  Entity  Data  Pile 
--  logical  urn  t  12. 


OPEN  ( UNI T=1 2 ) 

UR. I  fE  (  1 3 ,  1  20j ) 

1303  FORMAT  (  ’O’  jTL’S,  ’Measured  Drier  g  y  Data  File’,/, 


eO  100  1-1,1000 


Transfer  a  record  from  the  input 
file  to  the  Duffer  SCuRE. 


READ  (12, 1201 ,END= 101 )  uCORE 
1201  FORMAT  (A72) 

UR  I  f  E  (13,1 30 2* )  SCOPE 

C  -  - - if  SCORE (1:1)  Is  a  then  SCORE 

C  should  contain  data  and  a  valid 

keyword.  So,  see  it  SCORE <2: 4) 
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Program 


c: 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c. 

c 


e - - 


SUBROUTINE  PUTRND 


Veriicm:  84.11.30 

Author:  Paul  H.  Gstdiek 

Air  force  Institute  of  i  tchntiloyy 
Wr ight-Patterson  Air  Force  Base,  OH 

Descf  iption:  TFiis  routine  sloi  ts>  the  random  number  gener  ator  seed 

values  used  by  MINUM  i  ri  a  file  (unit  14)  for  future 
reference. 


SUBROUTINE  , UTRND  ilU,  IX) 
INTEGER-*  a  IU,  IX 


C  - Put  the  current  random  number 

C  generator  seed  values  into  a  tile 

C  for  future  reference. 


OPEN  (UNIT-14) 

WRITE  (14,1401) 
WRITE  (14,1402) 
WRITE  (14,1403) 
WRITE  (14,1401) 


WRITE 

(14, 1404) 

IU 

WR I TE 

(14, 140b) 

IX 

l  -i  0  1 

FORMAT 

( 7  r  • **********’  )  ,  ' 

1  **’  ) 

14c- 

FORMAT 

(  'if** 

Seed  s 

for  a  random  number  generator  in  M1NUM' , 

+  * 

***’  ) 

1403 

FORMAT 

(  ’  *** 

called  b/  program  DIATOM 

i  ’ 

P  *P '  ) 

1  404 

FORMAT 

(  *  > 1 U= ’  , 

115) 

1  40  a 

FORMAT 

(  ’  >!>•=’  , 

115) 

C  - - Write  the  end-of-file  mark  and 

C  c 1 ose  the  f i re. 


ENDFILE  (UN1T=14> 

CLOSE  (UNIT- 14, STATUS-’ KEEP’ ) 

RETURN 

END 
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I  >  em 


aUi ROUTINE  WnVL 


l'i 

c 

C  Vlt  ■_  ion:  34.11. 30 

c 

C  Author  :  Paul  H.  OstdieE 

Va 

0  Air  Force  Institute  oi  Technology 

C  Wr ight-Pattcrson  Air  Force  Euue,  OH 

C 

C  Description:  iinVE  control-  the  execution  ot  FUN  ta.no  therefore 

—  CI-eEr.)  ,  I  iuKi'inL  ,  and  PUT WAV  to  L%i  1  CU~ 

C  U'.t  iu'id  store  t  tie  nortna  1  i  zed  wave  functions.  This 

ie  done  in  an  iterative  fashion  for  each  wave  func¬ 
tion  once  the  eigenvectors  are  returned  from  FUN. 

C 


OUBkOUTINL  WAVE 
CM«kACTER*?2  eilbl 

iNIEGER  ELMT,  1,  NOEVAL,  30BN,  N,  NELMT,  EVMCNT,  EVMLVL(2b), 
♦  NOPNTS,  LEVEL,  HUNDRD 

klAL  E  VM  (  2b )  ,  STPS2E,  NODE  <  101),  PSIVAL(202>, 

+  BEGIN,  EVAL ( 202 )  ,  EVEC < 202 , 202 ) ,  H(2050u),  S<202,4), 


EVMWC5),  RES  I D  ( 2b )  ,  MIN,  HBAR,  MU,  L(202,4> 


C  OMI'loM 

/ENERGY/ 

EVH, 

EVMW ,  RESID, 

EVMLVL ,  EVMCNT 

COMMON 

/El GENS/ 

NELMT 

EVAL, 

,  BEGIN,  HBAR, 
EVEC,  NOEVAL, 

MU 

JOBN 

COMMON 

/MATRIX/ 

H,  S, 

L,  N 

common 

/  WAVFUI 1  / 

NODE, 

PSIVAL,  NOPNT 

S,  LEVEL 

C  ‘J  MM  ON 

/CHRLBL / 

3TLBL 

C  - - Call  function  FUN  using 

C  MI NUM’s  "best"  potential  energy 

C  parameter  set  to  get  the 

C  eigenvalues  and  vectors. 

MIN  -  FUN  (PAP i't) 

WRITE  (13,1301)  MIN 

'.301  f ORMAT  (’OThe  ^utn  ot  residuals  squared  is:  ',01b. 7) 

NODE  < 1 )  -  BEGIN 
I/O  10  1-2,  (NELMT+I  ) 

NODE ( I  )  =■  NODE ( I -I  )  +  STPS2E 

10  CONTINUE 


NOPNTS 

HUNDRD 


N 

100 


T1 


ELMT  =  25 

IF  (NGEVAL  .LT.  25)  ELMT  -  NGEVAL 

C  - Open  the  wave  function  output  file. 

OPEN  (UNIT  -17) 


C  - Write  the  1  aDe  1  /  cornmeri  t  read  in  from 

C  unit  12  in  IVEADER. 


WRITE  (17,1701)  STL EL ( i : 65) 

1701  FORMAT  ( * >LAEL= ' , AE5 ) 

C  - - - - - Write  Lhe  matrix  S  to  the  wave 

function  output  file. 


00  50  I  =  1  ,  N 

WRITE  (17,1702)  I,  (S(I,J),J=1,4) 

1^02  FORMAT  ( ’ >5 ’ , 15, 4615 . 7 ) 

50  CONTINUE 

C  - --Normalize  each  eigenvector  (wave 

C  function)  stored  in  EVEC .  Then 

C  write  it  to  the  wave  function 

C  output  file. 

DO  100  LEVEL= 1 , ELMT 
DO  110  1=1, N 

PSIVAL(I)  =  EVEC < I , LEVEL) 

110  CONTINUE 

C  --  - NORMAL  normalizes  the  eigenvectors 

C  (wave  function  values). 


CALL  NORMAL 


C  - PUTWAV  writes  the  wave  function  to 

C  the  output  file. 

CALL  PUTWAV 

100  CONTINUE 

C  - Write  an  End-Of-File  mark  arid  close 

C  the  output  file. 


EMDF1LE  (UNIT-17) 
CLOSE  ( UN  I T= 1 7 ) 

9999  RETURN 
END 


F-24 


n  non  non  n  o 


Pr  ogru.ii: 


SUBROUT I NE  NORMAL 


C 


I 


C 

c 

c 

c 

c 

c 

c 

c 


Version:  8-5.11.30 

Author:  Paul  H.  OstdieL 

Air  Farce  Institute  of  Technology 
Wright-Patterioii  Air  Fcrce  Best,  OH 

Description:  This  routine  normalizes  the?  wave  function  referenced 

by  LEVEL. 


SUBROUTINE  NORMAL 


INTEGER  I,  J,  K,  KK,  N,  NGEVmL,  NOPNTS,  LEVEL,  JOBN, 
EVMLVL<25),  EVMCNT ,  NELMT,  IE R,  II,  13,  1202 

REAL  EVM i 25)  ,  STPS2E ,  NODE < 101 ) ,  PSI VAL ( 202 )  , 

i  Lu  lu,  EVAH202),  LVEC  <202,202)  ,  H(20S03),  S<202,4), 
E  VHvJ  i  25 )  ,  RES  ID  (25)  ,  AREA,  SUM, 

A (202),  HEAR,  MU,  L (202,4) 


COMMON  /ENERGY/  EVM,  EVMW,  RESID,  EVMLVL ,  EVMCNT,  STPS2E, 
<  NELMT,  BEGIN,  HEAR,  MU 

COMMON  /El  GENS/  EVAL ,  EVEC ,  NOEVAL,  JOBN 
common  /matrix/  h,  s,  l,  m 

COMMON  / WAVFUN /  MODE,  PS  I VAL,  NOPNTS,  LEVEL 


Dr.TA  II,  13,  1202  /I,  3,  202/ 


AREA  =  0.0 


Multiply  SUPS I VAL  =  A 


VMULGF  is  an  IMSL  routine  for  matrix 
mu ltipl ication. 

C  (Band  Symmetric  storage  mode  times 

C  Full  storage  mode.) 

C Ai_.L  VMULQF  (S,  N,  13,  1202,  PSIVAL,  II,  1202,  A,  1202) 

c  - Multiply  PSIVAL*A  =  AREA 

C 

C  VMULFF  is  an  IMSL  routine  for  matrix 

C  multiplication. 

C  (Full  storage  mode  times  Full 

'•  C  storage  mode.) 


CALL  VMULFF  (PSIVAL,  A,  II,  N,  II,  II,  1202,  AREA, II,  1EK) 


c 

so 

9999 


Ncriiidlize  the  wave  function 


DO  50  I  =  1 ,  N 

PSIVAL(I)  =  PSIVAL(I)  /  SORT <Afc$< AREA) ) 
CONTINUE 

RETURN 

END 
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e - 

c 

C  Projram:  SUBROUTINE  PUTWAV 

t.. 

C  Version:  £*4.11.30 

C 

C  Author:  Paul  H.  Obtdiek 

C 

C  Air  Force  Institute  of  Technology 

C  Wr ight-Patterson  Air  Force  Ease,  OH 

C 

C  Description:  PUTWAV  writes  the  value  of  the  current  wave-function 

C  the  node  value,  and  the  cubic  spline  coef  -f  ic  lerits 

C  to  a  -file  (logical  unit  17). 

C 

c - 


SULROUTIME  PUTWAV 


INTEGER  NOPNTS,  LEVEL,  VIBLVL,  I,  J 


REAL  MODE (101),  PSIVAL(202> 

COMMON  / WAVFUN /  NODE,  PSIVAL,  NOPNTS,  LEVEL 

C  - Write  the  vibrational  quantum  number 

C  ,  a  count  number,  wave  function 

C  value  to  each  record. 

VIBLVL  =  LEVEL  -  1 


DO  10  I - 1 , NOPNTS 

WRITE  (17,1701)  VIBLVL,  I,  PSIVAL(I) 


l^Ol 

FORMAT 

1<-. 

CONTINUE 

vy  9  V 

RETURN 

END 
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Pr  ogr  din 


SUBROUT I WE  OUTPUT 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


Veraiori!  3*4.11.30 

Author:  Paul  H.  QstdieL 

Air  Force  Institute  of  Technology 
Wr i gh t -Patter son  Air  Force  Base,  OH 

Description:  OUTPUT  controls  the  execution  of  output  routines 

0ETP0T ,  PRNTER,  PLTRES,  and  PlTPOT. 


SUBROUTINE  OUTPUT  l PR ,  PP ) 

INTEGER  PR,  PP 

PEAL  R(IOOO),  POT VAL l 1000 ) 


C  - Subroutine  GETPOT  calculates  the 

C  value  of  the  potential  energy  model 

C  at  lOOO  grid  points. 


CALL  GETPOT  (R,  PGTVAL) 


C  - Subroutine  PRNTER  writes  to  the 

l  output  listing  file. 


CAl.L  RENTER  l  R ,  POT VAL) 


C  - - - Use  PLTRES  to  create  a  residual  plot 

C  file  if  the  user  set  PR=1  m  the 

C  input  file  (unit  11). 


IF  (PR  . ME .  0)  THEN 

CALL  PLTRES 
END  IF 


C  - - - Use  PLTPOT  to  create  a  potential 

C  energy  plot  file  if  the  user  set 

C  PP-1  in  the  input  file  (unit  11). 


IF  (PP  .ME.  0)  THEN 

CALL  PLTPOT  (R,  POTVAL) 
END  IF 


9999  RETURN 


C  Program:  SUBROUTINE  GET POT 

C 

C  Version:  8*1.11.30 

C 

C  Author:  Paul  H.  OstdieL 

C 

C  (iir  Force  Institute  of  Technology 

C  Wr ight-Pattersun  Air  Force  Base,  OH 

C 

C  Description:  This  routine  calculates  the  value  of  the  potential 

S  energy  model  at  1000  grid  points  using  the  best 

C  parameter  set. 

C 
C 


SUBROUTINE  GETPOT  (ft,  POTVAL) 

INTEGER  EVI-1CNT,  I,  MELMT ,  EVMLVL(25) 

PEAL  BEGIN,  STPS2E,  R(1000),  POTVAL  (  1 000 )  ,  RO,  Rl,  PO,  BP 0, 
*  PI,  DP  1  ,  EVM  ( 25 )  ,  EVNW  ( 25 )  ,  RESID(25>,  STEP,  HBA1-  MU 

COMMON  /ENERGY/  EVM,  EVMW,  RESID,  EVMLVL,  EVMCNT,  STPS2E, 

-  NELMT,  BEGIN,  HBAR,  MU 

STEP  =  (STPSL'E  *  NELMT)  /  1000.0 


DO  10  1-1,999 

RO  BEGIN  +  (STEP  *  (1-1)  ) 
Rl  =  RO  +  STEP 


C  - POTENT  returns  the  value  of  the 

C  potential  function  at  RO  and  Rl. 

CALL  POTENT  (RO,  Rl,  PO,  DPO,  PI,  DPI) 

R  (  I  )  -  RO 

POTVAL ( I )  =  PO 
IF  (I  . EQ .  999)  THEM 
R  <  I  +  1  )  =  Rl 

POTVAL  (IM)  =  PI 
END  IF 

10  CONTINUE 

9979  RETURN 
EMI) 
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Program:  SUBROUTINE  PRNTLK 

Version:  84. II. SO 

Author:  Paul  H.  Qstdiek 

Air  Force  Institute  of  Technology 
Ulr  i  gh  t -Pat  ter  son  Air  Force  Base,  OH 


DcbC i  i p  t i on : 


This  routine  prints  the  best  set  of  potential 
parameters,  the  value  of  the  potential  function 
using  these  parameters,  the  observed  and  calculated 
energy  levels  with  residuals. 


SUBROUTINE  PRNTER  (R,  PuTVAL) 

INTEGER  I,  J,  N,  NOEVAL,  PARl'lNO,  EVMLUL  ( 25 )  ,  EVMCNT, 
i-  NELMT,  IT1,  IT2,  IT3,  IT4,  ITS,  ITS,  STOP,  JOBN 

REAL  EVAL ( 202 ) ,  EVEC < 202 , 202 > ,  H(20503), 

►  S  (  202 , 4 )  ,  R(IOOO),  POTVAL ( 1000) ,  PARMUO),  PARMLO(IO), 

r  PARMH I (10) ,  EVM ( 25 ) ,  EUMW<25),  RESID(25), 

i-  STPS2E ,  BEGIN,  CONST(IO),  HEAR,  MU,  L(202,4) 

COMMON  /ENERGY/  EUM,  EVHW,  RESID,  EVMLVL,  EVMCNT,  STPSZE, 
i-  NELMT,  BEGIN,  HEAR,  MU 

COMMON  /FARMS/  FARM,  PARMLO,  PARMHI ,  CONST,  PARMNO 
COMMON  /E I GEMS/  FVAL,  EVEC ,  NOEVAL,  JOBN 
COMMON  /MATRIX/  H,  S,  L,  N 


-Rescale  the  parameters  for  display. 


DO  10  1=1, PARMNO 

PARM(I)  =  PARMLO ( I )  ♦  PARM ( I )  *  PARMHI (1) 

CONTINUE 


- Print  the  best  parameter  set. 

WRITE  (13,1301) 

FORMAT  ( ’ IThe  "best"  set  of  potential  energy  parameters  is:’) 

WRITE  (13,1302)  < PARM < I ) , I = 1 , 1 0 ) 

FORMAT  <’0’,5(G1S.7),/,5<G15.7)> 


-Print  the  value  of  the  potential 
function  at  1000  grid  points. 


WRITE  (13,1303) 

FORMAT  (’O’, ’The  values  of  the  potential  are:’) 
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DO  20  1=1,200 

WRITE  (13,1304)  R ( I ) , POTVAL ( I ) , R ( I +200 ) , POTVAL ( I +200 ) , 

<•  R(  1+400)  ,  POTVAL  (  I  +400)  ,R(  1+600)  ,  POTVAL  ( I  +  oG0 )  , 

+  R( 1+800) .POTVAL (1+800) 

1304  FORMAT  (’  ’ , 5 ( ’ > ’ , F? . 4 , 0 1 5 . 7 , 2X ) > 

20  CONTINUE 

C  - Print  the  observed  and  calculated 

C  energy  levels  and  the  residual. 

WRITE  (13,1309) 

1309  FORMAT  ('IThe  weighted  energies,  eigenvalues,  and  residuals:’) 

DO  25  1=1,25 


WRITE  (13,1308)  EVMLVL (  I  )  ,  EVMW ( I > ,  EVM ( I ) ,  EVAL  ( 1  )  ,  RESID(I) 
1308  FORMAT  (’  ’ , 12, 2X, F7. 5, 3G15. 7) 

25  CONTINUE 

C  - Print  the  calculated  energy  levels 


WRITE  (13,1305) 

1305  FORMAT  (’OThe  calculated  eigenvalues  are:’) 

STOP  =  24 

IF  (JOBN  .NE.  0)  STOP  =  N 

DO  40  1=0,24 
I  T 1  =1  +  1 

IT2  =  2a  +  I 
IT3  =51+1 
IT4  =76+1 
IT5  =  101  +  I 
IT6  =  126  +  I 

IF  (STOP  .LT.  IT6)  THEN 
IF  (STOP  .LT.  IT5)  THEN 
IF  (STOP  .LT.  IT4)  THEN 
IF  (STOP  .LT.  ITS)  THEM 
IF  (STOP  .LT.  IT2)  THEN 
IF  (STOP  . GE.IT1)  THEN 

WRITE  (13,130a)  I,  EVALII+1) 

FORMAT  (’  ’ , I3,2X,G15. 7> 

END  IF 
ELSE 

WRITE  (13,1307)  I,  EVAL < I + 1 ) ,  (1+26),  EVALd+27) 

FORMAT  ('  ’ ,2(I3,2X,G15. 7) ) 

ENDIF 
ELSE 

WRITE  (13,1311)  I,  EVAL ( I ♦ 1 ) ,  (1+26),  EVAL<I+27), 

(1+51),  EVAL ( I ♦ 52 ) 

FORMAT  (’  ’ ,3(I3,2X,G15.7) ) 

ENDIF 
ELSE 


1  3. -a 


1307 


1311 
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1310 

1312 

1312 

40 


9999 


WRITE  (13,1310)  I,  EVAL(I+1),  (I+2to),  EVALd+27), 

+  (I+Sl),  E0AL<l+52),  d+7o>,  EVrti.d+77) 

FORMAT  (’  •  ,4  d3,2X,G15.  7)  ) 

END  I F 
ELSE 

WRITE  (13,1312)  I,  KVALd  +  l),  d+26),  EVAL(I+27), 

■t  <  1  +  51 )  ,  E0AL(I+S2),  (  I  +  76 )  ,  EUAL(I+?7>, 

t  d  +  101),  EOAL  (  I  +  1 02 ) 

FORMAT  C  ’ , 5 (13, 2X, 015. 7) > 

END  IF 
ELSE 

WRITE  (13,1313)  I,  EUALd  +  1),  d+26),  EOALd+26), 

+  (1+51),  E0AL(I+52),  (1+76),  E0AL(I+77), 

+  d  +  101),  E0AL(I+102),  (1+126),  E0AL(I  +  127) 

FORMAT  (’  ’ ,6(13, 2X,G 15.7) ) 

END  IF 
CONTINUE 

IF  (STOP  .GE.  26)  THEN 
I  =•  25 

WRITE  (13,1306)  I,  EUAL(26) 

CNDIF 

RETURN 

END 
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■  V  f  V 
.  --  < 


' .  '  o  I,-  L-  . 


n  n 


c — - 

c 

C  Program:  SUBROUTINE  PLTRES 

C 

0  Version:  34. 11. SO 

C 

C  Author:  Paul  H.  Ostdiek 

C. 

C  Aii'  i-'orce  Institute  u+  Technology 

C  Wr l gh t -Pst terson  Air  Force  Base,  OH 

Description:  PLTRES  creates  a  -file  (logical  unit  15)  containing 

C  the  residuals  of  the  least  squares  Fit  for  plotting. 

C 


SUBROUTINE  PLTRES 

INTEGER  EVMCNT,  EVMLVL<25),  NELMT 


REAL  BEGIN,  EVM<25),  EVMW(25),  RESID125),  STPS2E,  HBAR,  MU 

COMMON  /ENERGY/  EVM,  EVMW,  RESID,  EVMLVL,  EVMCNT,  STPS2E , 

+  NELMT,  BEGIN,  HBAR,  MU 

C  - Open  the  output  file. 


OPEN  ( UN  l  T  =  1 5  ) 


C  - - - - Write  the  vibrational  quantum  number 

C  and  the  residual  value  to  the  file. 

DO  10  1=1, EVMCNT 

WRITE  (15,1501)  EVMLVL(I),  RESIDU) 

1501  FORMAT  (*  ’, 13,015.7) 

10  CONTINUE 

C  - Write  an  End-Of-File  mark  and  close 

C  the  output  file. 


ENDFILE  ( UNI T= 1 5 ) 
CLOSE  ( UNI T= 1 5 ) 

9999  RETURN 
END 
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C  Program:  SUBROUTINE  PL TROT 

C 

C  Version:  84.11.30 

C 

C  Author:  Paul  H.  Qstdiel; 

C 

C  Air  Force  Institute  ot  Technology 

C  Ur i ght  -Patter son  Air  Force  Base,  OH 

C 

C  Di.-3i:ript.Qn:  PLTPOT  creates  a  Pile  containing  the  value  of  the 

C  potential  energy  model  at  lOUO  grid  points. 

C 

c - 

SUBROUTINE  PLTPOT  (R,  POTVAL) 

INTEGER  I 

REAL  R(IOOO),  POTVAL (1000) 


C  - Open  the  output  Pile. 

OPEN  l  UN  I T  —  1  <s  ) 

C  - Write  the  grid  position  and  the 

C  value  of  the  potential  energy  model. 

00  10  1=1,1000 

WRITE  (lo, 1 P01)  R ( I ) ,  POTVAL ( I ) 

U,0  1  PORWAT  <’  ’,2015.7) 

10  CONTINUE 

- Write  an  End -Of -File  mark  and  close 

C  the  P i 1 e . 

ENDFILE  ( UN I T  = 1 6 ) 

1 1  CLOSE  ( UN I T= 1 6 ) 


Program:  SUBROUTINE  TKAlLk 

Vt-r  l.  i  on :  8A  .11. 30 

Author:  Paul  H.  UstdieK 

Air  Foret-  Inst  i  tuts  c+  Technology 
U right -Patterson  Air  Force  Base,  OH 

i)ticr  ip*,  ion:  TKAlLK  closes  the  output  listing  tile  (logical 

units  13  and  o )  and  stops  the  CPU  and  wall  time 
use  statistics. 


SUBROUTINE  TRAILR 


Shut  down  the  run  statistics. 


CALL  ETIME 
CALL  WTIME 


Close  the  listing  outputs. 


CLOSE  (UMIT= 13, STATUS=’ KEEP’ ) 
CLOSE  (UNIT =6) 


RETURN 


Wave-function  grids  do  not  agree’) 


C 

C 


9999 

*ADD 

*AOD 

*ADD 

TADD 

■-TADD 


ELi  E 

WRITE  (13,1302) 
FORMAT  COERRGR  -- 
ENDIF 


Close  the  output  Tile  and  run 
statistics. 


CALL  TRLFCF 
END 

HDRFCF 

RDRFCF 

CLCFCF 

OUTFCF 

TRLFCF 


n  n 


-Read  in  all  data  concerning  the 
wave  functions  of  the  upper  state 
from  logical  unit  12. 


LUNIT-12 

CALL  RDRFCF  (LUNIT) 


-Make  sure  the  grids  used  for  both 
states  have  the  same  number  of  grid 
points. 


IF  (NOPNTA  .  EC'.  NOPNTB)  THEN 


-Make  sure  Loth  wave  function  sets 
used  identical  S  matricies. 


TRUE . 


50 

1  =  1, 

NOPNTA 

IF 

(  SA  (  I 

,  1 )  . NE . 

SB  < I , 1 )  ) 

SOK  = 

. FALSE 

IF 

(SA  <  I 

,2)  .NE. 

SB< 1,2) ) 

SOK  = 

.FALSE 

IF 

(SA  <  I 

, 3 )  . NE . 

SE  < I , 3)  ) 

SOK  = 

.FALSE 

IF 

<  SA  (  1 

,  4  )  . NE . 

SB  < 1 , 4)  ) 

to 

o 

* 

II 

.FALSE 

IF 

<  .  NOT 

.  SOK)  GO 

i  TO  51 

CONTINUE 


-Calculate  the  Franck-Condon  factor 
for  every  combination  of  the  wave 
functions  from  the  two  states  and 
store  the  value  in  FACTOR (LVLA, LVLB) 


IF  (SOK)  THEN 
00  200  1=0,24 
LOLA  =  I 
DO  100  J=0,24 
LVLB  =  J 


--CLCFCF  computes  the  Franck-Condon 
factor  tor  the  two  wave  functions 
PSIA(LVLA)  and  PSIB(LVLB). 


CALL  CLCFCF  (LVLA,  LVLB) 


CONTINUE 

CONTINUE 


-0UTPCF  writes  a  table  of  Franck 
-Condon  factors  to  the  listing. 


CALL  0UTFCF 


ELSE 

WRITE  <13,1301) 
FORMAT  < ’ 0ERR0R 
ENDIF 


-  Matricies  S  do  not  match’) 
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Appendix  H 
Program  FCFACT 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


Program:  FCFACT 

Version:  84.11.20 

Author:  Paul  H.  Ostdiek 

Air  Force  Institute  of  Technology 
Wr lght-Patterson  Air  Force  Base,  OH 

Description:  This  program  calculates  the  square  of  the  inner 

product  between  each  wave  function  from  two  sets. 
This  is  the  Franck -Condon  factor.  Each  wave 
function  set  is  expected  to  have  25  (v=0  to  24) 
wave  functions  in  the  format  used  by  DIATOM. 

I/O  logical  unit  6  --  output  listing  file 

11  --  input  file  (wave  set  1  v’> 

12  --  input  file  (wave  set  2  v“) 

13  --  output  listing  file 

(same  as  unit  6) 


PROGRAM  MAIN 
CHARACTER*72  LBLA,  LLLB 

INTEGER  I,  J,  LULA,  LVLB,  LUNIT,  NOPNTA,  NOPNTB 
LOGICAL  SOK 

REAL  FACTOR ( 0 : 24 , 0 : 24 ) ,  PG I A < 1 : 202 , 0 : 24 >  ,  PS I B ( 1 : 202 , 0 : 24 ) , 

+  SA ( 202 , 4 ) ,  SB ( 202 , 4 ) 

COMMON  /DTAFCF/  FACTOR,  PSIA,  PSIB,  SA,  SB,  NOPNTA,  NOPNTB 
COMMON  /LABELS/  LBLA,  LBLB 


Open  the  output  listing  file,  print 
the  header,  and  start  run  stats. 


CALL  HDPFCF 


C  Read  in  all  data  concerning  the 

C  wave  functions  of  the  lower  state 

C  from  logical  unit  11. 


LUNIT  *  11 

CALL  RDRFCF  (LUNIT) 


H-l 


Appendix  G 
Program  FCFACT  Flow 


G-l 


FALSE 


FIRST  = 
ENDIF 


PREFIX  =  CONST ( 1 ) 
RATIO  CONST  (  2) 


(TP (2)  -  TP ( 1 ) )  * 
Ri 


-Calculate  the  value  o+  the 
potential  energy  PI  and  its  slope 
DPI  at  the  right  edge  Rl  of  the 
grid  el emen t . 


TP (3) f (PREFIX* (TP  <21* (RAT IG**TP( 1 ) ) -TP ( 1 ) * ( RATI 0**TP ( 2 1 ) )  ) 
PREFIX*TP(1)*TP<2> / (TP (2) -TP ( 1 ) ) * (RATIO**TP ( 2 > -RAT IO**TP ( I ) ) 


-Limit  infinity  to  a  usable  number 


IF  (PI  .GT.  1E10)  THEN 
PI  =  IE  10 
DPI  = - 1 E 1 0 
ENDIF 


-Calculate  the  value  of  the 
potential  energy  PO  and  its  slope 
DPO  at  the  left  grid  edge  RO.  Also 
limit  infinity. 


IF  (PI 

PO  ■ 
DPO  ■ 


10000.0) 


ENDIF 

ELSE 

RATIO  =  CONST (2)  /  RO 

PO  =  TP (3) + (PREFIX* (TP ( 2 ) * ( RAT I 0**TP ( 1 ) ) -TP ( 1 ) * ( RAT 10* * TP ( 2 ) ) ) ) 
DPO  =  PREF I X*TP <1 ) *TP (2) / (TP (2) -TP ( 1 ) ) * (RaTIO**TP (2) -RATI0**TP ( I ) ) 

IF  (PO  .GT.  1E10)  THEN 
PO  =  IE  10 
DPO  --1E10 
ENDIF 
ENDIF 


RE  .’URN 


n  r> 


C  program:  SUBROUTINE  POTENT  (Mle) 

Version:  84.11.30 

C 

C  Author:  Paul  H.  Gstdieh 

C 

C  Air  Force  Institute  of  Technology 

C  Ur i ght -Pat ter  son  Air  Force  Base,  OH 

C 

C  Description:  This  routine  returns  the  potential  energy  value  and 

C  slope  at  the  left  and  right  grid  element  boundaries. 

C  The  Mie  function  is  used  a  model. 

C 

C - 

SUBROUTINE  POTENT  (RG,  Rl,  PO,  DPQ,  PI,  DPI) 

INTEGER  1,  PARMNO 
LOGICAL  FIRST 

REAL  PARl'l  (10),  PARMLO  (10),  PARMHI  <  10)  ,  RO,  Rl  ,  PO,  DPO,  PI,  DPI, 

+  TP ( 1 0 )  ,  CONST  1 1 0) ,  PREFIX,  RATIO,  UPPER110) 

COMMON  /FARMS/  FARM,  PARMLO,  PARMHI,  CONST,  PARMNO 

DATA  FIRST  / . TRUE. / 


C  - Scale  the  potential  parameters  PARM 

C  to  get  correct  parameter  values  TP. 


DO  10  1=1, PARMNO 

TP  (  I  )  =  PARI’ILO  (  I  )  +  PARM(I)  *  PARMHI  (  I  ) 

UPPER ( I )  =  PARMLO! I)  +  PARMHI  II) 

CONTINUE 
IF  (FIRST)  THEN 
URITE  (13,1301) 

FORMAT  (/, /,  ’OPoter.tial  Model  used  is  Mie’> 

WRITE  (13,1302) 

FORMAT  (’02  Constants  arid  3  Parameters  are  used’,/, 

+  ’  Constant  1  is  the  Dissociation  Energy’,/, 

r  ’  Constant  2  is  the  I n ter nuc 1  ear  Separation’,/, 

♦  ’  Parameter  1  is  the  power  alpha’,/, 

+  ’  Parameter  2  is  the  power  beta’,/, 

*  ’  Parameter  3  is  the  energy  shift’) 

WRITE  (13,1303) 

FORMAT  (’0  Number  Constant  Parameter  L’, 

i  ’ower  Limit  Upper  Limit’,/,’  -  -’, 

,  *  ’  ,3  (  ’  ’)) 

DO  20  1=1,10 

WRITE  (13,1304)  I,  CONST(I),  TP ( I ) ,  PARMLO i I ) ,  UPPER ( 1 ) 

FORMAT  (’  ’ ,5X, 1 2 , 5X , 4 ( 2X , G 1 b . 7 > ) 

CONTINUE 

F-4l 


10 

lG'Jl 

1302 

i  :.*(  3 

1  304 
20 


END  IF 


PREFIX  CONST  (1)  *  TP  (2)  i  (TPil)  -  TP  (2)) 

RATIO  -  CONST (2)  /  R1 

IF  (RATIO  .GT.  6.5)  RATIO  -  o.5 

- Calculate  the  value  of  the 

potential  energy  PI  and  its  slope 
DPI  at  the  right  edge  R1  of  the 
grid  element. 

PI  =  TP (3) + (PREFIX* (RATI G**TP ( 1 ) - ( (TP ( 1 ) /TP (2) >*RA1 IO**TP (2)) ) ) 
DPI  -  PREFIX  *  ( TP ( 1 > /Ri)*( RAT  1 0**TP i 2 )  -  RAT  I Q**TP ( i )  ) 


Limit  infinity  to  a  usable  number. 


IF  (PI  .GT.  1 E 1 0 )  THEN 
PI  =  1E10 
DPI  = - 1 E 1 0 
END  IF 


Calculate  the  value  Df  the 
potential  energy  PO  and  its  slope 
DPO  at  the  left  grid  edge  RO.  Also 
limit  infinity. 


IF  (RO  .EG.  0)  THEN 

IF  (PI  .GE.  10000.0)  THEN 
PO  =  PI 
DPO  =  DPI 
ELSE 

PO  =  1E10 
DPO  =- 1 E 1 0 
END  IF 
ELSE 

RATIO  -  CONST (2)  /  RO 

IF  (RATIO  .GT.  6.5)  RATIO  =  6.5 

PO  =  TP (3) + (PREFIX* (RAT I 0**TP ( 1 ) - ( ( T P ( 1 ) / TP ( 2 ) ) *RA T IO**TP ( 2 ) ) ) ) 
DPO  =  PREFIX* ( TP ( 1 ) /RO ) * ( RATI G**TP ( 2 ) -RATIQ**TP ( 1 ) ) 

IF  (PO  .GT.  IE 10)  THEN 
PO  =  1 E 1 0 
DPO  = - 1 E 1 0 
END  IF 
END  IF 


RETURN 


Pr  ogr  am : 


SUBROUTINE  POTENT  (Lennar d-Jonei) 


Version:  84.11.30 

Author:  Paul  H.  Qstdiek 

Air  Force  Institute  o-f  Technology 
■fir  lght -Patterson  Air  Force  Base,  OH 

Description:  This  routine  returns  the  potential  energy  value  and 

slope  at  the  le+t  and  right  grid  element  boundaries 
The  Lennar d -Jones  model  is  used. 


SUBROUTINE  POTENT  IRO,  Rl,  PO ,  DPO,  PI,  DPI) 

INTEGER  I,  PARMNO 
LOGICAL  FIRST 

REAL  HARM  (10),  PARHLO(IO),  PARMHK10),  RO,  Rl,  PO,  DPO,  PI,  DPI, 
TP ( 1 0 )  ,  C  ONST (10),  PREFIX,  RATIO,  UPPER(IO) 

COMMON  /FARMS/  FARM,  PARMLO,  PARMHI ,  CONST,  PARMNO 

DATA  FIRST  / . TRUE. / 


Scale  the  potential  parameters  FARM 
to  get  correct  parameter  values  TP. 


DO  10  1=1, PARMNO 

TP ( I )  -  PARMLO ( I )  ♦  FARM ( I )  *  PARMHI (I) 

UPPER ( I )  =  PARMLO ( I )  +  PARMHI  (I) 


CONTINUE 
IF  (FIRST)  THEN 
WRITE  (13,1301) 

FORMAT  (/,/,’ OPotent  ial  Model  used  is  Lenard -Juries '  ) 
WRITE  (13,1302) 

FORMAT  (’02  Constants  and  3  Parameters  are  used’,/, 

’  Constant  1  is  the  Dissociation  Energy ’ , / , 


’  Const u  n  t 

2  l  s 

the 

1  n ter nuc 1 ear  Separation’,/, 

’  Parameter 

1  1  * 

the 

power  a 1 pha ’ , / , 

’  Pat  a  me  ter 

2  i  * 

the 

power  beta ’  , /  , 

’  Pc *  r  Sine  t  e  r 

2  l  & 

the 

energy  shift’ ) 

WRITE  (13,1303) 

FORMAT  (’0  Number  Constant  Parameter  L’, 

’ awer  Limit  Upper  Limit’,/,’  -  -’, 

. - '  ,  3  (  ’  - ')) 

DO  20  1=1,10 

WRITE  (13,1304)  I,  CONST (I),  TP ( I ) ,  PARMLO ( 1 )  ,  UPPER  l  1 ) 

FORMAT  (’  ’ ,bX, I2,bX,4(2X,GlS.7) ) 

CONTINUE 


FALSE 


FIRST  =  . 

END  IF 

EXPNET  -  TP (  1  )  *  (CONST (2)  -  Ri ) 

PREFIX  =  2  *  CONST ( 1 )  *  TP ( 1 )  *  EXP(EXPNET) 


- /-Calculate  the  value  of  the 

potential  energy  PI  and  its  slope 
DPI  at  the  right  edge  Ri  of  the 
grid  element. 

PI  -  CONS  I  (1)*(  (  <  1-EXP  (EXF'MET)  1**2)  -I)  *TP(2> 

DPI  =  PREFIX* ( 1-EXP (EXPNET) ) 


Limit  infinity  to  a  usable  number 


IF  (PI  .or.  1E10)  THEN 
PI  1E1G 
DPI  =-lL.O 
END  IF 


Calculate  the  value  of  the 
potential  energy  PO  and  its  slope 
DPO  at  the  left  grid  edge  RO.  A1 
limit  infinity. 


IF  (RO  . EQ.  0)  THEN 

IF  (PI  .GE.  10000.0)  THEN 
PO  =  PI 
DFO  =  DPI 
ELSE 

PO  -  1E1G 
DPO  =~ 1 E 1 0 
ENDIF 
EL  CE 

CXPNET  =  TP { 1 )  *  (CONST (2)  -  RO) 

PREFIX  =  2  *  CONST ( 1 )  *  TP ( 1 )  *  EXP(EXPNET) 

PO  -  CONST ( 1 >*(< (1-EXP (EXPNET) ) **2) -1 > +TP (2) 
DPO  =  PREFIX*( 1 -EXP (EXPNET) ) 

IF  (PO  .GT.  1E10)  THEN 
PO  =  1E10 
DPO  =  -  IE  1 0 
ENDIF 
ENDIF 


RETURN 


c-  - 

c 

C  Program:  SUBROUTINE  POTENT  (Morse) 

C 

C  Version:  84.11.30 

C 

C  Author:  Paul  H.  Ostdiek 

C 

C  Air  Force  Institute  of  Techno  1 agy 

C  Wr  lght-Pattersori  Air  Force  Base,  DH 

C 

C  Dcscript ion:  This  routine  returns  the  potential  energy  value  and 

C  slope  at  the  left  arid  right  grid  element  boundaries. 

C  The  Morse  function  model  is  used. 


SUBROUTINE  POTENT  <R0,  Rl,  PO,  DPO,  PI,  DPI) 
INTEGER  I,  PARMNO 
LOGICAL  FIRST 


REAL  FARM (10),  PARMLO(IO),  PARMHI ( 10) ,  RO,  Rl,  PO,  DPO,  PI,  DPI, 
+  TP (10),  CONST (10),  PREFIX,  EXPNET,  UPPER (10) 

COMMON  /PARMS/  PARM,  PARMLO ,  PARMHI ,  CONST,  PARMNO 

DATA  FIRST  /.TRUE./ 


C  - Scale  the  potential  parameters  PARM 

C  to  get  correct  parameter  values  TP. 


DO  10  1=1, PARMNO 

TP ( I )  =  PARMLO ( I )  +  PARM ( I )  *  PARMHI (I) 

UPPER ( I )  =  PARMLO ( I )  +  PARMHI (I) 
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1301 

1302 


1303 


1304 

30 


CONTINUE 
IF  (FIRST)  THEN 
WRITE  (13,1301) 

FORMAT  </,/,’ OPotent l a  1  Model  used  is  Morse’) 

WRITE  (13,1302) 

FORMAT  (’02  Constants  and  2  Parameters  are  used’,/, 

♦  ’  Constant  1  is  the  Dissociation  Energy’,/, 

+  ’  Constant  2  is  the  In ternuc 1  ear  Separation’,/, 

+  ’  Parameter  1  is  the  factor  beta’,/, 

«■  ’  Parameter  2  is  the  energy  shift’) 


WRITE  (13,1303) 

FORMAT  (’0  Number  Constant  Parameter 

♦  ’ qhci  Limit  Upper  Limit',/,’  - 

.  . - ’  ,  3  (  ’  - ’)) 

DO  20  1=1,10 


WRITE  (13,1304)  I,  CONST ( I ) ,  TP ( I ) ,  PARMLO ( I ) ,  UPPER < I ) 
FORMAT  <’  ’ ,5X, I 2 , 5X , 4 ( 2X , G1 5 . 7 ) ) 

CONTINUE 


F-39 


I-’  , 


- Calc U late  the  value  of  the 

potential  energy  PI  arid  its  slope 
DPI  at  the  right  edge  R1  of  the 
grid  element. 

PI  =  T P  < 1 )  X  Pi** 2  i  TP  (2) 

DPI  -  2  *  TP ( 1 )  *  R 1 

- - Limit  infinity  to  a  usable  number. 

IF  (PI  .GT.  1E10)  THEN 
Pi  =  1E10 
llPl  -  -  1 E  1 0 

END  I F 

- - Calculate  the  value  of  the 

potential  energy  PO  and  its  slope 
DPO  at  the  left  grid  edge  RC.  Also 
limit  infinity. 

PO  =  TP<1)  *  R0**2  TP  (2) 

DPO  =  2  *  TP ( 1 )  *  RO 

RETURN 


C  Program:  SUBROUTINE  POTENT  (Single  Harmonic  Gsc i 1 1 dtor ) 

C 

C  Version:  34. 11. SO 

C 

C  Author:  Paul  H.  Ostdiek 

C 

C  Air  Force  Institute  of  Technology 

C  Wr i gh t -Patterson  Air  Force  Ease,  OH 

C 

C  Descr ipt ion:  This  routine  returns  the  potential  energy  value  and 

C  slope  at  the  left  and  right  grid  element  boundaries. 

C  The  Harmonic  Oscillator  model  is  used. 

C 

c - 

SUBROUTINE  POTENT  (RO,  Rl,  PO ,  DPO,  Pi,  DPI) 

INTEGER  I,  PARMNO 
LOGICAL  FIRST 

PEAL  PARM (10),  PARHLG (10),  PARMHI (10),  RO,  Rl,  PO,  DPO,  PI,  DPI, 

♦  TP (10),  CONST (10),  PREFIX,  EXPNET,  UPPER l 10) 

COMMON  /FARMS/  PARM,  PARMLO,  PARMHI ,  CONST,  PARMNO 

SAT  A  FIRST  / . TRUE.  / 


C  - - Scale  the  potential  parameters  PARM 

C  to  get  correct  parameter  values  TP. 


DO  10  1=1, PARMNO 

TP ( I )  =  PARMLG(I)  +  PARM ( I )  *  PARMHI ( I ) 

UPPER! I)  =  PARMLO ( I )  +  PARMHI ( I ) 

CONTINUE 
IF  (FIRST)  THEN 
WRITE  (13,1301) 

FORMAT  </,/,’ ©Potential  Model  used  is  Harmonic  Oscillator’) 
WRITE  (13,1302) 

FORMAT  ('00  Constants  and  2  Parameters  are  used’,/, 

+  ’  Parameter  1  is  the  power  alpha’,/, 

+  ’  Parameter  2  is  the  energy  shift’) 

WRITE  (13,1303) 

FORMAT  (’0  Number  Constant  Parameter  L’, 

+  ’ower  Limit  Upper  Limit’,/,’  -  -’, 

♦  . - .  >3(  »  - *  )  ) 

DO  20  1=1,10 

WRITE  (13,1304)  I,  CONST ( I )  ,  TP ( I >  ,  PARMLO ( I ) ,  UPPER  (  I ) 

FORMAT  (’  ’  ,5X, I 2 , SX , 4 ( 2X , 0 1 S . 7 )  ) 

CONTINUE 
FIRST  =  .FALSE. 

END  IF 
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Only  one  of  the  remaining  four  subroutines  is  used  as 
subroutine  POTENT.  The  following  are  routines  for  the  Single 
Harmonic  Oscillator,  Morse,  Lennard- Jones ,  and  Mie  potential 


r 


C 

C  Program:  SUBROUTINE  HDRFCF 

C 

C  Version:  84.ii.S0 

c 

C  Author:  Paul  H.  OstdieL 

C 

C  Air  Force  Institute  of  Technology 

C  Wr lght-Patterson  Air  Farce  Base,  OH 

C 

C  Description:  HDRFCF  opens  the  output  listing  Tile  (logical  units 

C  13  and  6)  and  starts  CPU  and  nail  time  use 

C  statistics. 

C 

C - 

SUBROUTINE  HDRFCF 

CHARACTER*3  VERSN 
CHARACTER# 13  PCN 

INTEGER*3  IDATE (3) ,  ITIMEC3),  IUSER ( 4 ) 

VERSN  =  ’ 84. li . 30* 

PCN  =  ’ GEP/84D-6/ 1 . 4 ’ 

C  - Initiate  the  run  statistics. 

CALL  BT I HE 
CALL  ST  I ME 


C  - Get  the  current  date,  time,  and 

C  user  name  Tor  output  on  the  header 


CALL  DATE (IDATE) 
CALL  TIME ( IT  I ME) 
CALL  USERNOt IUSER) 


Open  the  output  listing  file  and 
write  out  the  header 


OPEN  (UNIT-13) 

OPEN  ( UNIT-6 ) 

WriTE  (13,1301)  IUSER,  VERSN,  IDATE,  PCN,  I T I  ME 
1  FORMAT  (’1  User:  ' , 4A3 , T5 1 ,  ’Air  Force  Institute  of  Technology’, 

•  T 1 1 0 , ’Version:  ’ , A8 , 

•  /,’  Date:  ’ ,3A3,T114, ’PCN:  ’,«13, 

/,’  Tune!  ’  , 3h3, T57,  ’  FRANCK-CONDON  FACTORS’,/,/) 
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Program:  SUBROUTINE  RDRFCF 

Version:  84.11. 30 

Author:  Paul  H.  Ostdiek 

Air  Force  Institute  of  Technology 
Wr i ght-Patter son  Air  Farce  Base,  OH 


Desc  r i p  1 1  on : 


This  routine  reads  all  records  From  one  oF  two  input 
Files  (LUNIT-il  or  12).  Each  record  read  is  written 
to  the  output  listing  (unit  13).  Data  records  are 
marked  by  a  ’ > ’  in  column  1.  These  records  contain 
data  referenced  by  a  single  key  word.  All  other 
records  are  considered  comments.  This  routine  uses 
an  internal  read,  eg  READ  (SCORE (5: 19) , ’ (E15. 7) ’ )  X 
reads  From  columns  S  to  19  of  the  character 
variable  SCORE  using  the  edit  descriptor  E15.7  into 
the  real  variable  X. 


SUBROUTINE  RDRFCF  (LUNIT) 

CHARACTER*?2  LBLA,  LBLB 
CHARACTER*80  SCORE 

INTEOER  I,  J,  LUNIT,  NOPNTA,  NOPNTB,  COUNT,  TSTLVL ,  SCNT 

REAL  FACTOR (0: 24, 0: 24) ,  PS I A ( 1 : 202 , 0 : 24 ) ,  PS I B ( 1 : 202 , 0 : 24 ) , 
►  SA ( 202 , 4 ) ,  SB (202, 4) 

COMMON  /DTAFCF/  FACTOR,  PSIA,  PSIB,  SA,  SB,  NOPNTA,  NOPNTB 
COMMON  /LABELS/  LBLA,  LBLB 


-Open  logical  unit  LUNIT. 


OF  EM  (UNIT-LUNIT) 


DO  100  1-1,10000 


-Transfer  a  record  From  the  input 
File  (unit  11  or  12)  to  the  buffer 
SCORE. 


READ  (LUNIT, 1101, END-101)  SCORE 
FORMAT  (A80) 


-IF  SCORE (1:1)  is  a  *>’  then  SCORE 
may  contain  data. 
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IF  (3C0RE(1:1>  .  EQ.  ’ > ’ ) THEN 

IF  ( SCORE  (2:3)  .EG.  ’  LABL '  >  THEN 


This  record  contains  a  label  'for  the 
state  involved. 


IF 

(LUNIT  .Ed. 

11) 

LBLA ( 1 

:  72) 

=  SC0RE(7: 

78) 

IF 

(LUNIT  . £Q . 

12) 

LBLB ( 1 

:  72) 

=  SC0RE(7: 

78) 

ELSE 

IF 

(SCORE (2: 2) 

.  £0 . 

•  S’  ) 

THEN 

- This  record  contains  data  for  the 

£  matrix. 

READ  (SCORE (3:5) , 1 (13) * )  SCNT 
IF  (LUNIT  .EG.  11)  THEN 

REAL  (SCORE (6: 20) , ’ (E15. 7) * )  SA(SCNT,1) 

READ  ( SCORE  (21:35)  , ’  (E15.7)  ’  )  SAiSCNT,2) 

READ  (SCORE (36: SO) , ’ (£15.7) ’ )  SA(SCNT,3) 

READ  (SCORE (51 : 65) , ’ (E15. 7) ’ )  SA(SCNT,4) 

ELSE 

READ  (SCORE (6: 20) , ’ (E15. 7) * )  SBtSCNT.l) 

READ  (SCORE (21 : 35) , ’ (E15. 7) ’ )  SB(SCNT,2) 

READ  (SCORE (36:50) (E15. 7) * )  SB(SCNTf3) 

READ  (SCORE (51 : 65) , ’ (E15. 7) * )  SB(SCNT,4> 

ENDIF 

ELSE 


This  record  must  have  data  For  the 
TSTLVL  th  wave  function. 


READ  (SCORE (2: 5) , ’ ( 14) * )  TSTLVL 
READ  ( SCORE (o:9),'(I4)’)  COUNT 
IF  (LUNIT  .£0.  11)  THEN 

READ  (SCORE ( id: 24) , ’ (E15. 7) ' ) 
IF  (COUNT  .GT.  NOPNTA)  NOPNT  A 
ELSE 

READ  (SC0KE( 10: 24) , ’ (E15. 7) ' ) 
IF  (COUNT  .GT.  NGPNTB)  NOPNTB 
END  IF 
END  IF 
END  IF 
END  IF 
CONTINUE 


PS I A ( COUNT , T  STLVL  > 
=  COUNT 

PS I B (COUNT , TSTLVL ) 
=  COUNT 


Close  logical  unit  LUNIT. 


CLOSE  (UNIT  =  LUN I T ) 


RETURN 


C  Program:  SUBROUTINE  CLCFCF 

C 

C  Version:  8*4.1.1.30 

C 

C  Author:  Paul  H.  Ostdiek 

C 

C  Air  Force  Institute  of  Technology 

C  Wr 1 gh t -Patter son  Air  Force  Base,  OH 

C 

C  Description:  This  routine  calculates  the  Franck-Condon  factor  for 

C  the  ^ave  functions  referenced  by  LOLA  and  LOLB. 

C 

c - 

SUBROUTINE  CLCFCF  (LOLA,  L0l_B) 

INTEGER  I,  J,  LUNIT,  NOPNTA,  NOPNTB,  LOLA,  LOLB,  IER, 

*  II,  13,  1202 

REAL  FACTOR  (0:  2*4 , 0 :  2*4  )  ,  PSI A  !  1 :  202 , 0 :  24 >  ,  PSI B  (  1 :  202 , 0 :  2*4  )  , 

^  SA(202,4),  SB  < 202 ,  *4 )  ,  A(202),  0A1202),  0B(202), 

r  AREA 

COMMON  / DT  AFCF /  FACTOR,  PSIA,  PSIB,  SA,  SB,  NOPNTA,  NOPNTB 

DATA  II,  13,  1202  /l,  3,  202/ 

DO  100  I -1, NOPNTA 

OA ( I )  =  PS  I A  <  I , LOLA ) 

OB ( 1 1  =  PSIB (I, LOLB) 

100  CONTINUE 

AREA  =0.0 


C  - Multiply  SA*0A  =  A 

C 

C  OMULGF  is  an  IMSL  routine  for  matrix 

C  multiplication. 

C  (Band  Symmetric  storage  mode  times 

C  Full  storage  mode.) 

CALL  OMULQF  ( SA ,  NOPNTA,  13,  1202,  OA,  II,  1202,  A,  1202) 

C  - Multiply  0B*A  =  AREA 

C 

C  OMULFF  is  an  IMSL  routine  tor  matrix 

C  multiplication. 

C  (Full  storage  mode  times  Full 

C  storage  mode.) 


CALL  VMULFF  <VB,  A,  II,  NOPNTA,  II,  II,  1202,  AREA,  II,  IER) 


c  - The  Franck-Condon  factor  is  AREA 

C  squoued. 

FACTOR  (LVLA.LVLB)  =  ARLA**2 

9999  RETURN 
EMD 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


Program:  SUBROUTINE  OUTFCF 

Vers  i  dm  :  84.li.30 

Author:  Paul  H.  Ostdiek 

Air  Force  Institute  of  Technology 
Wr ight-Patterson  Air  Force  Base,  OH 

Description:  This  routine  prints  a  25 

table  (to  unit  13). 


by 


25  Franck-Condon  factor 


SUBROUTINE  OUTFCF 
CHARACTER*72  LBLA ,  LBLB 
INTEGER  I,  J,  NOPNTA,  NOPNTfi 

REAL  FACTOR (0:24,0:24) ,  PS I A ( 1 : 202 , 0: 24 ) ,  PS IB ( 1 : 202, 0: 24 ) , 

*  SA(202,4),  SB  (202,4) 

COMMON  / DTAFCF /  FACTOR,  PSIA,  PSIB,  SA,  SB,  NOPNTA,  NOPNTB 
COMMON  /LABELS/  LBLA,  LBLB 

WRITE  (13,1305)  LBLA,  LBLB 

1305  FORMAT  (’  v "  (across  the  page  -  lower  state)  is  for  ’,/, 

+  ’  1 , A72 ,/,/,’  v*  (down  the  page  -  upper  state)  is  for’,/, 

+  ’  ’ , A72 ) 


WRITE  (13,1301)  ( J , J -0 ,15) 

1301  FORMAT  ( ’ 0 ’  ,  * v  ' \ v  “  '  , 12,  1 5 ( 5X ,12),/) 

C  - Write  the  first  part  of  FCF  table. 

DO  10  1=0,24 

WRITE  (13,1302)  I,  ( FACTOR ( J , I ) , J =0 , 1 5 ) 

1302  FORMAT  (’  ’,I2,1X,16(2X,F5.3)) 

10  CONTINUE 


WRITE  (13,1303)  (J,J=16,24> 

1303  FORMAT  ( *  1 ’  , 1 X , 9 ( 5X , 12 )  , / ) 

C  - Write  the  second  part  of  FCF  table. 

DO  20  1=0,24 

WRITE  (13,1304)  I,  ( FACTOR ( J , I ) , J = 1 6 , 24 ) 

1304  FORMAT  (’  ’ , 12, IX , 9 ( 2X, F5. 3 ) ) 

20  CONTINUE 

9990  RETURN 


ill 


n  n 


C  Program:  SUBROUTINE  TRLFCF 

C 

C  Version:  8A. 11.30 

Author:  Paul  H.  Qstdiel; 

C 

C  Air  Force  Institute  of  Technology 

C  Wr i gh t -Patterson  Air  Force  Base,  OH 

C 

C  Description:  TRLFCF  closes  the  output  listing  file  (logical 

C  units  13  and  6)  and  stops  the  CPU  and  wall  time 

C  use  statistics. 

C 

- - 

SUBROUTINE  TRLFCF 

C  - Shut  down  the  run  statistics. 

CALL  ET I ME 
CALL  WTIME 

C  - Close  the  listing  outputs. 

CLOSE  ( UN IT— 1 3 , STATUS^ ' KEEP ’ ) 

CLOSE  ( UNI T=6 ) 


999*?  RETURN 


Program  Flow  Symbols 


A  process  of  some  kind. 


A  module,  e.g.  SUB1  may 
be  a  subroutine  or  a 
function. 


YF3 

>  A  decision  point. 


) 


Start  or  stop  a  task. 


A  disk  file. 


Output  listing  (printer) 


Off  page  connector 


Appendix  J 


The  Single  Harmonic  Qsci 1 lator 

The  Schrodinqer  wave  equation  describing  a  one 
dimensional,  single  harmonic  oscillator  is  (10:75): 


+  %E  -  |kx2)*  =  0  (J-l) 

dx-  TT 


where 


|kx2  =  i[l(02x2  (J-2) 

Bounded  solutions  exist  only  for  the  discrete  energy  levels 
defined  by  the  vibrational  quantum  number  v  and: 

E(V)  =  ftw(v  +  5)  (J-3) 

This  equation  shows  that  these  energy  levels  are  spaced 
equally  apart. 

The  orthonormal  wave  functions  of  the  harmonic 

oscillator  are  then  described  in  terms  of  Hermite  polynomials 
1 

Hv(a2x)  as: 


where  the  normalization  constant  Nv  is 


and  the  scaling  constant  a  is 

a  =  1AH 
a  -h 

The  first  ten  Hermite  polynomials  are: 

H0(q)  ■  1 
Hl(q)  =  2(1 
H2(q)  =  4<>2  '  2 

H3(q)  =  8<13  -  12« 

H4(q)  =  16q<*  '  4842  +  12 

H5(q)  =  32q5  "  l60q3  +  120q 

H6(  j  =  64q6  -  480q4  +  720q2  -  120 

j  =  I28q7  -  1344q5  +  3360q3  -  l680q 

H8(q)  =  256<18  ‘  3584q6  +  l3440q4  -  13440q2  +  1680 

H »  =  51 2q9  -  92l6q7  +  48384q3  -  80640q3  +  30240q 
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Paul  H.  Ostdiek  was  born  in  Ardmore,  Oklahoma  on  13 
June  1957*  He  graduated  1'rom  Park  Hills  High  School,  Fair¬ 
born,  Ohio  in  1975-  He  received  appointments  to  the  U.S.  Air 
Force  Academy,  U.S.  Naval  Academy,  and  U.S.  Coast  Guard 
Academy  that  same  year.  In  1979  he  graduated  from  the  U.S. 
Air  Force  Academy,  and  was  commissioned  a  second  lieutenant 
in  May  of  that  year.  His  first  assignment  was  to  the  princi¬ 
pal  laboratory  of  the  Air  Force  Technical  Applications  Center 
at  McClellan  AFB.  There  he  was  a  Scientific  Data  Systems 
Analyst  his  first  year,  and  Chief,  Resources,  Computer 
Operations  his  second.  During  his  third  year  he  was  Chief 
of  the  Microanalysis  section.  During  his  fourth  year  he 
designed  and  implemented  a  Scientific/Management  Data  Base 
System  for  the  laboratory.  In  November  1981  he  received 
the  Air  Force  Commendation  medal.  He  was  selected  Company 
Grade  Officer  of  the  Quarter  for  McClellan  AFB  in  June,  1982. 
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Institute  of  Technology. 
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This  thesis  developed  a  finite  element  solution  of  the 
Schrodinger  wave  equation.  This  technique  is  used  by  a 
computer  program  to  calculate  the  energy  levels  and  wave 
functions  of  a  diatomic  molecule  for  a  particular  poten¬ 
tial  energy  model.  The  potential  energy  model  is  a  function 
cf  a  set  of  parameters  which  a  non-linear  minimization  rou¬ 
tine  varies  before  solving  the  wave  equation.  This  is  done 
in  an  iterative  manner  until  the  calculated  energy  levels 
agree  in  a  least  squares  sense  with  the  observed  energy 
levels.  Then  the  transition  probabilities  (Franck-Condon 
factors)  between  the  wave  functions  are  calculated  by 
another  program  developed  for  this  thesis.  Finally,  two 
programs  were  written  to  determine  the  energy  levels 
observed  in  spectroscopic  data.  One  uses  Dunham  coefficients 
and  the'-^Dunham  equation  while  the  second  uses  a  least  square 
fit  to  the  data  directly. 

The  four  programs  were  tested  and  appear  to  work  correct- 
Ly.  The  numeric  solutions  were  compared  with  the  analytic 
solutions  of  the  single  harmonic  oscillator.  The  lowest 
25  energy  levels  agreed  to  within  0.005$  accuracy  while 
their  wave  functions  appear  to  agree  to  within  0.40 $  i 

accuracy. 
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